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Cos3flaH B lesItxX HHPOPMUPOBaHHA YUTATeIbCKOM ayAMTOPMU O HOBeMMWHX JOCTWKEHHAX HM MepcneKTHBax B OOACTH 
MeXaHHKH, MalIMHOCTpoeHHA, MHPOpMaTHKH MW BEIMMCIMTeIbHOM TexHuKH. V3yqaHue saBiaetca opyMoM JUIA 
COTpyHHYecTBa POCCHMCKHX HU MHOCTpaHHBIX YYCHBIX, CIOCOOcTBYyeT COMKEHHIO POCCHicKOrO H MHpoBoro Hay4Ho- 
MHOpMallHOHHOrLO MpocTpaHcTBa. 


2KyphHasl BKUOUeH B Tepedenb pelleH3HpyeMbIX HAYUHbIX 131aHhi, B KOTOPOM JOJDKHbI ObITh ONyOJIMKOBAaHbI 
OCHOBHbIe Hay4Hble pe3y.IbTaTbl AUcceprauHi Ha COMCKaHHe y4ueHon cTeneHu KaHHAaTa HayK, Ha cOonCcKaHHe 
yueHor cTeleHu OKTOpa Hayk (lepeyens BAK) no ceaxyrouuM Hay4HbiM CciewMasIbHOCTAM: 


1.1.7 — Teoperwueckas MexaHika, JMHaMHkKa Maly (TexHH4eCKHe HayKH) 

1.1.8 — Mexannka Jedopmupyemoro TBepyoro Tesla (TeXHH4ecKHe, (PH3HKO-MaTeMaTHYecKHe HayKH) 

1.1.9 — Mexanuka 2kHJKOCTH, ra3a H Wla3MbI (TexHu4eCcKHe HayKH) 

1.2.2 — MarematTw4eckoe MojesIMpoBaHHe, YMCICHHbIC METOBI H KOMIVICKChI IporpaMM (TexHHYecKHe HayKH) 

2.3.1 — CucTemubiit anamu3, ynpaBsenve u oOpaboTKa HHopMallMH, CTaTHCTHKa (TexHM4eCKHe HayKH) 

2.3.3 — ABTOMaTH3al[HA HM yIpaBeHHe TeXHOJIOrM4eCKHMH TIpoleccaMH HM MpOH3BOACTBaMH (TeXHM4ecKHe HayKH) 

2.3.5 — Matematw4eckoe 1 IporpaMMHoe OGeciieveHHe BLTMHCIHTeJIBHBIX CHCTEM, KOMIUICKCOB H KOMMBIOTEPHBIX ceTeli (TeXHH4eCKHe HayKH) 
2.3.7 — KommbiotepHoe MoyesMpoBaHve VM aBTOMAaTH3al[HA IpOeKTHpOBaHHA (TeXHHYeCKHe, (PU3HKO-MaTeMaTHYecKHe HayKH) 

2.3.8 — Uacdopmatuka 4 HHC*OpMallHOHHEIe Mpoleccsl (TexHu4eckve HayKH) 


2.5.2 — Mammmuopeyenue (TeXHHYeCKHe HayKM) 

2.5.3 — Tpeuue vf 43HOC B MalliMHax (TeXHM4eCKHe HayKH) 

2.5.5 — Texnonorna 4 o6opyqoBaHve MexaHHyeckoli H PU3MKO-TexHHYecKol OOpaboTKH (TexHH4eckHe HaykKH) 

2.5.6 — TexHosIOrMa MalMHOCTpoeHHA (TeXHHYeCcKHe HayKH) 

2.5.8 — CBapka, poCTBeHHBbIe MIpoleccbl U TeXHOMOrMN (TexHHyeckHe HayKH) 

2.5.9 — Metogpl 4 MpHOopbI KOHTpOsIA MW WMarHOCTHKM MaTepuHasion, H3]evIMi, BeELeCTB H IpHpogHO cpebI (TexHM4ecKHe HayKH) 
2.5.10 — Pugpapmuueckne MalIMubl, BaKyyMHat, KOMMpeccopHasd TeXHuKa, THApO- H WHeEBMOCHCTeMBI (TEXHHYeCKHe HayKH) 
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Abstract 


Introduction. The speed rise of railway transport and an increase in the loads on the axles of wheelsets necessitate the 
modernization of the existing fleet. Scientific studies in the field of rolling stock dynamics are aimed at taking into account 
the oscillatory processes that occur during the movement of railway vehicles in a traditional design. The attachment of 
supplementary elements was considered at the coupling level of two cars and the attachment of a third trolley in the center 
of gravity of the railway platform. The scientific literature has not paid enough attention to the construction of 
mathematical models that make it possible to assess the dynamic states of such constructive solutions. The objective of 
this study is to create a method for evaluating the dynamic conditions of a car. The situation is considered when an 
additional set of mass-inertial and elastic elements is introduced into its structure, and the general dynamic condition of 
the vehicle depends on the adjustment of their parameters. 

Materials and Methods. The basic research tool is the structural mathematical modeling, which is based on an approach 
where the source design scheme is a mechanical oscillatory system in the form of a solid body on elastic supports with 
supplementary typical elements introduced into its structure. The dynamic analogue of the calculation scheme used is the 
block diagram of the automatic control system, the use of which provides detailing the connections between typical elastic 
and mass-inertia elements. 

Results. A method for estimating the dynamic states of railway vehicles is proposed. It is based on the construction of 
mathematical models, taking into account the introduction of an additional structure of mass-inertia and elastic elements. 
The impact of additional parameters on the dynamic condition of the vehicle is investigated. Analytical relations have 
been obtained that provide reducing the dynamic loads on the major structural elastic elements when changing the 
corresponding parameters of a technical object. The transfer function of interpartial relations is given, which provides 
controlling the interaction between the coordinates of the vehicle movement under the action of two kinematic 
disturbances of the in-phase type. 

Discussion and Conclusion. The generated mathematical model provides for assessment, monitoring and control of the 
dynamic state of the vehicle under conditions of kinematic disturbances. The research results can be used to modernize 


existing vehicles and create new ones with improved dynamics. 
Keywords: vehicle dynamics, dynamic condition, structural mathematical modeling, block diagram 
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AHHOTanHa 

Beedenue. Yaenuuenne ckopoctei TBWKeHHA 2%Kee3HOAOPOXKHOTO TpaHciopTa MW MOBbIMeHHe Harpy30K Ha OcH 
KOJICCHBIX Tap OOYcuaBIMBaIOT HEOOXOAMMOCTb MOJepHH3alHH CyujecTByromjero mapKa. HayyHble uccreqoBaHHA 
B OOaCTH JHHAMHKH MOABYWKHOTO COCTaBa HalipaBJIeHbl Ha yYEeT KOeOaTeMbHBIX MPOWeCcoB, BOZSHHKAIOWIMX IIpu 
JBWOKCHHM %Kee3HOLOPOXKHBIX TPaHCHOPTHBIX CpeACTB B TpayMUMOHHOM KOHCTpyKTHBHOM HcnomHeHuH. IIpuco- 
C2 MHEHHe JOMOUHUTCIBHBIX 3JICEMCHTOB PaCCMaTPHBAJIOCh Ha YPOBHE CIeIKH JBYX BarOHOB H IIpHcoeqMHeHHH Tpe- 
Tbe TEMOKKU B WEHTpe THKECTH %*Kee3HOMOpOXKHON maTpopMpt. [locrpoeHuto MaTeMaTHYeCKHX MOJeeH, 103BO- 
JIAIOUIMX OLWCHUTh WHHAMMYeCKHe COCTOAHHA TAKHX KOHCTPYKTHBHBIX pelieHHi, B Hay4HOH WMTepatype He yeueHoO 
ocTaTOuHO BHUMaHHA. Llenb JaHHOrO HCceOBAaHHuA — CO3aTb METOA OWCHKH AMHAMH4eCKHX COCTOAHHH BaroHa. 
PaccmaTpuBaeTca CHTYala, KOra B ErO CTPyKTypy BBOAMTCA JOMOIHUTeIbHaA COBOKYMHOCTb Macc-HHepUMOHHBIX 
MW YIpyrux 9ICMCHTOB, IIpH¥eM OT KOPpeKTHPOBKM UX MapaMeTPOB 3aBHCHT OOUIee WHHaMMYecKOe COCTOAHHE 
TPaHCIOPTHOTO CcpeCTBa. 

Mamepuaavi u memoooi. ba30BbIM HHCTPyMCHTOM IIpoBeACHHA HCCIEOBaHHH ABJIAeTCA CTPYKTYpHOe MaTeMaTH- 
YeCKOe MOJCIMPOBAaHHe, B OCHOBE KOTOPOLO JIeKUT NOAXOA, Kora UCKOAHAaA pacueTHad CXeMa IIpecTaBIAeT COOOH 
MeXaHHyeckKylo KOJIeOaTeIbHy!IO CHCTeMy B Bue TBEPAOrO Tella Ha yIpyrux onopax C AOMONHUTeIbHOM BBe]EeH- 
HBIMH B e€ CTPyKTYpy THMOBbIMH 3IeMeHTaMH. J[MHAaMH4eCKHM aHaIOrOM HCIOb3YeMOM pacueTHOHM CXeMBI ABJIA- 
eTCA CTpyKTYpHasd CXeMa CHCTeMbI AaBTOMATHYCCKOTO ypaBICHHA, IPHMeHeEHHe KOTOPOM MO3BOsIAeT WeTaIM3upo- 
BaTb CBA3H Me2KAY THIOBLIMH YIIPyrHMH MW MacC-HHepIHOHHBIMH 3JIEMeHTaMH. 

Pe3zyiavmamot ucciedosanua. WpeynoxKeH MeTOA OWCHKM AHHAMMYCCKHX COCTOAHHM 2KeIIe3HOMOPOXKHBIX TpaHc- 
TIOPTHBIX CpeJCTB, OCHOBaHHBIM Ha MOCTPOCHHH MaTeMaTHYeCKHX MOJeeH, C YHeTOM BBECHHA JOMOMHUTeIbHOM 
CTpyKTypbI Macc-HHepIMOHHBIX HW yupyrux semMeHToB. UccneqoBaHo BIMAHHE JOMNONHUTCIBHBIX TapaMeTpoB Ha 
JMHaMHYeCKOe COCTOAHHE TpaHcnopTHoro cpeycTBa. Tlosy4eHbl aHasyIMTHYeCKNe COOTHOMEHHA, MO3BOIAIOMNe pu 
VM3MCHCHHH COOTBETCTBYIONHX HapaMeTpoB TEXHHYeCKOTO OOBEKTA CHH3UTb JMHAMMYeCKHe Harpy3KH Ha OCHOBHBIC 
KOHCTPyKTHBHbIe yipyrnve 91eMeHTHI. IIpuBpeyena nepeqaTouHad PYHKUNA MexXKMapUMaIbHbIX CBA3ZCH, MO3BOIIAIO- 
Wad KOHTPOJMpOBaTbh B3aHMOeCHCTBHe MEXKAY KOOPAMHATaMHM ABMXKCHHA TpaHCHOpTHOro cpeACTBa pH AelicTBUU 
TByX KHHeMaTH4eCKHX BOSMyIeHH CHHa3Horo TMA. 

O6cyotcoenue u 3akniouenue. CcpopMupoBaHHad MaTeMaTHYeCKad MOJCJIb MO3ZBOACT OLCHHTh JHHAMUYECKOE CO- 
CTOAHHE 2Kee3HOAOPOXHOTO TpaHCHOpTHOTO CpeACTBa B YCNOBHAX eMCTBUA KMHEMATHYeCKHX BO3MyMeHuH. Pe- 
3YJIbTATbI HCCICOBaHHM MOTYT ObITb IpHMeHEHE! IPH MOAepHH3alHu CYLWIECTBYIOIWIMX M CO3{aHHH HOBBIX TpaHc- 


TIOPTHBIX CpeCTB C yalyuleHHori WMHaMUKOH. 


Ksro"ueBble C10Ba: AWHAaMHKa TpaHCiopTHbBlx CpeACTB, THHAMUYCCKOC COCTOAHHE, CTPYKTYPHOe MaTeCMaTHYeCKOe 


MOJeIMpOBaHne, CTpyKTypHaa cxema 


BaarogapHoctu. ABTopbl BbIpaxaloT OarofapHOCTb 3acilyKeHHOMYy JleaTemio HayKu PO, 7.T.H., mpodeccopy 
Ennceesy C.B., a Take peqakIHH UM pelleH3eHTaM 3a BHUMAaTeJIbHOe OTHOMEHNHE K CTaTbe HW YKa3aHHble 3aMedaHHA, 


KOTOpBbIle MO3BOJIMJIN MOBbBICHTb Ce KayecTBo. 
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Introduction. The expansion of technical-and-economic ties at the interregional level, providing the growth of the 
country's industrial engineering potential, and maintaining the growing system of international trade and commercial 
relations largely depends and relies on rail transport [1]. The rail traffic volume is constantly increasing. This drives 
the need to take into account the perverse effects of increased dynamic loads [2], which directly affect the reliability 
of operation of both rolling stock and the track structure. Despite such negative factors, it is necessary to fulfill the 
planned indicators, which include an increase in the service speed, compliance with weight and length standards for 
trains, increasing axle load to 30 tons or more. This reflects the real demands of the development of the Russian 
economy and stimulates the creation of new more powerful locomotives, the renewal of the fleet of rolling stock, and 
the modernization of track facilities [3]. At the same time, the possibility of negative consequences of the 
intensification of transportation processes should also be considered. One of the major issues is the increase in the rate 
of wear of the track structure with the corresponding resulting difficulties [4]. 

Currently, particular attention is being paid to the development of a methodology for assessing the dynamic 
condition of rolling stock, the interaction of technical means and rail tracks, energy savings, and improving the 
reliability and safety of transportation processes. The mathematical modeling methodology is described, e.g., in [5]. 
At the same time, there are other possibilities for finding rational solutions [6]. It is important to pay attention to the 
modernization of the existing fleet of freight wagons, whose operation is no longer effective under increased loads [7]. 
One of the approaches that could be adopted for the development is the concept of installing an additional two-axle 
bogie for freight 4ax wagons [8]. In this case, we can expect a more uniform distribution of the load on the track 
structure, as well as the possibility of increasing the weight of the transported goods while maintaining the regulations 
for axial load values within 22 tons [9]. The unevenness of the wheel-rail contact parameters initiates the oscillatory 
movements of the car, which in turn forms the oscillatory movements of the car body. The oscillation process is also 
formed by the conditions of interaction of wagons inside the train [10]. In this case, dynamic restraint forces occur, 
superimposing on the static components of the overall reaction, which can significantly increase the level of dynamic 
interactions in the wheel-rail contact [11, 12]. However, the possibilities of structural mathematical modeling in 
assessing the dynamic states of railway vehicles with the introduction of additional links have not yet been given 
enough attention. Therefore, the objective of this study is to form a method for evaluating the dynamic states of a 
railway vehicle when introducing an additional set of mass-inertial and elastic elements into its structure, the correction 
of the parameters of which would affect the overall dynamic state of the vehicle. 

Materials and Methods. The methodological basis of the research is the structural theory of vibration isolation 
systems, which provides for the accurate assessment of the dynamic properties of a railway vehicle in a linear 
formulation, taking into account concentrated parameters and small oscillations relative to the position of static 
equilibrium or a steady-state process. The design scheme is a mechanical vibratory system with dynamic equivalent in 
the form of a block diagram of an automatic control system. This makes it possible to detail the connections between 
the elements of the vehicle, as well as to use methods characteristic of the theory of automatic control (transfer 
functions, transformations of structural circuits, convolution and simplification, frequency characteristics) [13]. 

Standard design schemes of freight railway vehicles are described by well-known calculation schemes, and their 
dynamic features can be estimated using linear calculation schemes. The article considers a railway vehicle in the form 
of a four-axle freight wagon designed for the transportation of heavy goods such as coal, ore, sand, rolled metal, 


small-sized metal structures, etc. (Fig. 1). 
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4 4 /// 1/7 


Fig. 1. Schematic diagram of a four-axle freight wagon 


The structure of the presented railway vehicle contains a body with mass M and moment of inertia J. It is based on 
two four-wheel bogies, conventionally represented as a set of mass-inertial and elastic elements. The dynamic features of 
the system under consideration show the excessive impact on the structural elements of the bogies of the railway vehicle. 
In this regard, to improve its dynamic properties, a third two-axle truck is additionally introduced into the suspension 


structure, which is located in the center of the wagon (Fig. 2). 


J 
_d 


ame 


Fig. 2. Schematic diagram of a six-axle freight car 


With a car's own weight of about 23 tons, and an axial load of 22 tons, this can provide an increase in the weight 
of transported goods by approximately 20 tons, i.e., significantly increase the efficiency of transportation processes 
without creating excessive dynamic loads on the track structure (TS). The practical implementation of the proposed 
approach requires the modernization of the design of a standard wagon. The modernization consists in creating an 
additional bogie attachment unit to reduce the load on the axle of a freight wagon by introducing an additional two- 
axle track. The upgrading of the attachment unit is carried out in the same structural and technical forms as the fastening 
using pins in two “standard” two-axle bogies, i.e., through the installation of over- and underpin bolsters with the 
corresponding standard-type pin assembly. 

The proposed method of increasing the efficiency of using four-axle freight cars under conditions of additional 
loads is aimed at solving the problem of increasing the load on the axle of the wheelset to 30 tons, and speeding up the 
trains. This is achieved through upgrading a typical freight car by installing an additional two-axle bogie with an 
appropriate device that provides conditions for its dynamic interaction with the frame structure of the freight car. 

The installation of an additional two-axle bogie by redistributing the load between a common set of wheelsets 
provides the possibility of transporting heavy loads while reducing the load on the axle. This maintains better operating 
conditions for the track and the superstructure while maintaining an acceptable length of the train. 

It is assumed that the pin assembly has an elastic rubber gasket that provides cushioning for the interaction of the 
bogie and the frame of the wagon body. At the same time, it is not supposed to upgrade the wheels of the wheelsets. 
Only the form of wiring of the elements of the “regular” pneumatic braking system and the configuration of the supply 


of pneumatic pipelines are changed. 
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Research Results. The design scheme of the considered railway vehicle in the first approximation can be 
represented as a mechanical oscillatory system consisting of a solid body with mass m and moment of inertia J, based 


on three elastic elements with stiffness k1, ko, k2. Kinematic effects are represented by in-phase harmonic functions of 


the same frequency (Fig. 3). 


Yo M,J 


Fig. 3. Design scheme of a railway vehicle 


under kinematic disturbance (z1, z2, Zo) 


The center of mass of the system — point O is located at distances /; and /, from the ends of the solid (points A1, B1). 
The element with stiffness ko is fixed at points D and Dj; distance OD, is indicated as /9. The movement of the system is 


considered in coordinates y1, y2 and yo, @, associated with a fixed basis. The calculations use the following ratios: 


Yo =a +by,,b=c(y, VWiyeN = Vo “41, V2 = Vo tld, 
l, l, 1 (1) 


: c= : 
1, +1, [+l L+L 


Vp =Vo thd, a= 


Lagrange formalism is used to derive differential equations of motion [14], which requires the construction of 


expressions for kinetic and potential energies. In this case, we have: 


T= m(ay\ +by3)° +Je*(ys-yi)’, (2) 


1 1 1 
M=Sh (91-21) +>40(¥0 -20) +5 ha (¥2-22) (3) 


Note that potential energy (3) can also be written as: 


1 1 1 
II me (1 ~z,)’ +S (¥2 =z) oa [ay + by, +1yc(y> ~y)-29] = 


1 2 1 2 1 2 4) 
su (1 —2,) #78 (2 —Z2) 7 [ay +Diy2 —zo| 5 
where a; = a—loc, b} = b + he. 
The system of equations of motion in coordinates y, y2 in the time domain takes the form: 
VW (ma? 4 Je”) + a (x, + ka? )-y3 ( fe? —mab) + Vokayb, = kz, +koayzZo, (5) 
V5 (mb? + Je?) + y (ks +kob? ) -y, ide? —mab) + Vykoayby =kyZ. +k ob Zo. (6) 


After applying the Laplace integral transformations under zero initial conditions [15], the system of equations (5), (6) 
can be represented in the operator form: 


yy [ (ma? + Je?) p? +k, + koa |-F [ (4? —mab) p? — koayby | = k,Z, +koa,Z, (7) 
J, [ (mb? + Je?) p? +k + kyb? |- 7, [ (4? —mab) p? — kyayby | = kyZ, tkybiZp, (8) 


where p = j@ — complex variable (7 = V—1 ), indicator <—> above the variable means its Laplace image [5]. 
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On the basis of (7), (8), in accordance with the provisions of the method of structural mathematical modeling [5], a 
structural mathematical model in the form of a block diagram is constructed, which is dynamically equivalent to an 


automatic control system (Fig. 4). 


(Jc? — Mab) p* — k,a,b, 


1 1 
2 2 

_M = > 
(Ma? + Jc’) p* +k, +kga; We aH) NathP, (Mb? + Jc”) p* + ky + kody ie 

Z, Yi Z 

k, < + | k, le 2 

29 z, 
ka <—_— = ha 


Fig. 4. Structural mathematical model of a railway vehicle 


Discussion and Conclusion. The feature of the system is that it has two partial blocks, each of which determines the 


corresponding partial frequencies: 


k, + koa? 

2 Lt Koa) 

= ——| > 9 

ie ma? + Jc? 0) 
2 

oe ky + kobj (10) 


° mb? + Jc?” 
which, in turn, determine the boundaries of the location of the natural oscillation frequencies of the system as a whole: 
Dicos <i <NZ <O5e06, (11) 
where @1c06, ®2cos — natural oscillation frequencies of the system. When operating on them, resonant modes may occur. 


Among the features of the system there are several simultaneously acting external disturbances. Assuming, for the 


sake of simplification, that Z; =Z, =Z) (quite acceptable at the stages of preliminary dynamic estimates), we consider 
that force factors act on the input of the first and second partial blocks: 

O, =7(k, +koa,), (12) 

QO, =7(k, + kod, ). (13) 

Using the block diagram in Figure 4, we write down the expressions for the transfer functions, assuming that there is 


a relationship between the external perturbation factors formed by the ratio: 


0, =0, QO, =a, (14) 
W(p)=2= (ke, +koay)| (mb? + Je?) p? + ky + kob? |+0(ky + kobs)| (Je? mab) p? ~ hoa] a 
Z A(p) 
Ws(p)= Ya 7 a(k, + kb) | (ma? +Jc?) p? +k, + kya? | +(k, + kya;)| (Je? —mab) p? hoards] iss 
Z A(p) 
where 
A(p)= [ (ma? +Je?) p? +k, + kya? || (mb? +Je?) p? +k, + kop? | -[(ve? ~mab) p? — koayby | (17) 


is the frequency characteristic equation of the system. 
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The numerator of the transfer functions in expressions (15), (16) determines the modes of dynamic vibration damping, 


which can be detailed from the equations obtained by “zeroing” numerators (15), (16): 


ae (A, + kos )( ea + ko? )— a1 (ey + obs /koandy figs 
me (ky + kay )(mb? + Je? ) + o1 (ky + kobi)(Je? —mab)’ 


ee (Ks + Koby )( ki + koa? )— a1 (ky + koa) oad, 6) 
a (kik + hob, )(ma? + Je?) + (kyky + koa )( Je? — mab) 


It follows from expressions (18), (19) that in a system with two degrees of freedom, in the presence of two 
interconnected disturbing factors, dynamic vibration damping modes at two frequencies may occur, whose parameters 
depend on the values of connectivity coefficient a. This coefficient can take negative, positive and zero values. 

The analysis of the structural mathematical model (Fig. 4) also implies the possibility of a special dynamic mode at 
the frequency: 

kya,b 
one = 3p (20) 
when the partial blocks get the possibility of disconnection. In this case, the system (Fig. 3) splits into two autonomous 
blocks, which do not create situations of interaction of partial structures. 

The implementation of this mode can cause a significant difference in the values of deviations at points A; and B), and 
a “spread” of values of the dynamic reactions at points 41, B, and D,. The determination of dynamic reactions at points 
A, B and D can be implemented according to the methodology described in [14], in which the dynamic reaction is defined 
as the product of a dynamic displacement (e.g., points 41, B; and D)) by the value of the reduced dynamic stiffness. 


For coordinate y,, the dynamic offset is determined from the expression for transfer function (15), and for 
coordinate y, — from expression (16). The frequency characteristic equation is used to determine the reduced 
dynamic stiffness (17): 


(4? —mab) p? —kyayby | 


kya*a, 


Kwi(p) =A, + _ : (21) 
pi (P)= hs (mb? + Jc?) p? +k, +kobp (mb? +Jc?) p? +k, + kobe 
Similarly, the value of the reduced dynamic stiffness in coordinate yz can be found: 
2 
; | ; kob2b, [ (ve? —mab) p? — kyayby | 22) 
Ty Pp = + = . 
ep ' (ma? + Jc?) p? +k, +koa? (ma? + Jc?) p? +k, + kya? 


To determine the specific values Kept ( p) and Kenpo ( P); it is necessary to find the value of frequency w*, which is 


determined from the condition that y, / y, =1. More generally, it is assumed that: 
22 ay, (23) 


where y — coefficient of connectivity of the movement of elements by coordinates y, and y,. Thus, this coefficient 
(for the case y=1) can be written as: 
- a(k, + koby)| (ma? + Je) p? +k, +kga? ] +(k,+ koa,)| (Je? - mab) p? — koayby | 


= (24) 
Wn (k, + koa,)| (mb? + Je?) p? +k + kob? |+a(k, +koby)| (Je? —mab) p? — koayby | 


Taking specific values y, it is possible to find the oscillation frequencies for the motion of the object in question in 


coordinates y, and y,. Specifically, at y = 1, the frequency of translational vertical vibrations of a solid is determined 


from the expression: 
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Fy _ ake + koi (dy + koa? )— (ki + koa kab; + (ki + koa) koarby - 
yi (ky +k ob, )(ma? + Je?) +(k; + kya, )( Je — mab) - 

—(K, + kya, )(ko + kobe) — (hea + Kobi )(dr + koa?) + 01(ky + kod: Koad, 
—(k, + koa ) (mb? + Jc? ) —a(ky +kob )( Je? - mab) 


2 


@ 
(25) 


Frequency @ of in-phase harmonic disturbances z;(f), z2(¢) and zo(t) provides the motion of a solid body at » = 0, 


L.e., ¥j =V2 =p. For other values y (y 4 1), y,; and y, are determined, and on their basis, for geometric reasons, 
value Vp is determined. For a solid body with known y, and y,, the position of the center of rotation (or oscillations) 


can be easily found [12]. 


The dynamic constraint reactions R Ale R pz and R p1 can be found in the first approximation by the formulas: 
Ra =hy,, Ray =k,yo, Rp =kYp. 
In general, dynamic reactions are determined using dynamic displacements y,, y., Vp, that are defined by 
expressions (15), (16). As for the dynamic displacement at point D, expression yp =a,y,+b,¥2 is used, whose 


parameters are determined by the above values. In expressions (15), (16), data on the connectivity of force factors of 
influence (parameter a) can be entered. To obtain specific data on the values of dynamic reactions, frequency 


parameters are introduced at which the ratio of oscillation amplitudes y, and y, is implemented through the 


coefficient of coupling of oscillation amplitudes y. 

The overall restraint force at points A, B and D is determined by the sum of two components: static and dynamic. The 
static component can be found from the expression for the transfer functions of dynamic reactions at p = 0 and setting the 
parameters of the static load (weight of the wagon and the cargo being transported). With the intensive development of 


oscillatory processes, when fluctuations in coordinates ¥,, ¥7, Vp increase, the overall reaction can vary significantly and 


differ from the static restraint force. In the presence of a dynamic component, the overall reaction can take on various 
values, in particular, zero or negative. 

The mathematical model formed within the framework of the proposed method is indicated by expression (25). It 
provides assessing the dynamic state of railway vehicles when additional connections are introduced into their structure 
to form a set of recommendations for obtaining stable modes of operation. The investigation of the features of the system 
using approaches typical of structural mathematical modeling allows us to consider in detail the connections between the 
elements. With regard to the technical object in question, this makes it possible to adjust the dynamic state of the technical 
object, based on varying the parameters of the set of additionally introduced elements, to reduce the load on the main 
parts of the suspension, as well as to establish the presence of natural frequencies and frequencies of dynamic vibration 
damping in the system. 

In the future, it is planned to conduct research on the introduction of dampers and motion conversion devices into the 
structure of the vehicle to assess the possibilities of structural mathematical modeling. Another interesting area is the 
assessment of the possibilities of changing dynamic reactions depending on external actions, which will allow us to assess 


the efforts exerted on various elements of the vehicle suspension. 
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Abstract 


Introduction. Slender structures of subsea energy production systems are under constant influence of currents and waves. 
Hydrodynamic loads result from the interaction of subsea pipelines, umbilicals, equipment supports with fluid flows, and 
lead to the vortex formation in the area behind the structures. Vortex-induced forces are the sources of the cyclic loading. 
They accelerate gradually the fatigue damage, which may result in a failure. One of the ways to reduce the loads on subsea 
structures is to alter the shape of a cross-section, taking into account the flow regime. Dependence of the resulting 
hydrodynamic loads on the cross-sectional shape and relative position of structures has not been studied in details for the 
uniform flow in the critical mode. The current work is aimed at filling this gap. The research objective is to consider the 
impact of the distance between the structures, and also, the presence of a D-shaped structure, placed upstream relative to 
the group of three cylinders of different cross-sectional shapes. 

Materials and Methods. The computational fluid dynamics approach was used in this work for numerical simulations of 
vortex-induced forces in the ANSYS Fluent software for cylinder with D = 0.3 m. Modelling was conducted with the 
Detached Eddy Simulation (DES) method, which combined advantages of the Reynolds-averaged Navier-Stokes equation 
(RANS) method and the Large Eddy Simulation (LES) method. The object of the research was the system of four 
structures in the 2D computational domain, which included the upstream D-shaped cylinder and the main group of three 
cylinders with the circular, squared and diamond shapes of the cross-section. The transient process was considered, where 
structures were under the influence of the uniform flow in the critical regime at Re = 2.5x10°. 

Results. Five sets of data were obtained in simulations for the time-dependent coefficients of the lift and drag forces: for 
the main system — of the D-shaped, circular, square and diamond structures, and also for the four systems — of only 
D-shaped, only circular, only square and only diamond shaped structures. Additional analysis was conducted for the effect 
of the distance between the structures on the amplitude of fluctuating hydrodynamic force coefficients. The obtained 
results are presented as time histories of coefficients of the lift and drag forces, frequency analysis and contours of 
velocity, pressure and vorticity fields. The results indicate a positive effect of the upstream D-shaped structure on reducing 
the drag force, acting on the central structure in the group of three cylinders located downstream. 

Discussion and Conclusion. The results of the performed studies facilitate the informed decisions regarding the 
arrangement of subsea structures in a group of four objects, depending on the cross-sectional shape and the distance 
between the structures. The upstream D-shaped structure provides reducing the hydrodynamic drag force acting on the 
central structure in the downstream group of three structures, thereby slowing the fatigue accumulation and increasing 
the time of safe operation. 


Keywords: vortex-induced forces, drag coefficient, lift coefficient, uniform flow 
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AH&aJIM3 BO3MO2KHOCTH CHHDKEHHA JIOOOBOTO CONPOTHBJICHHA 3a cuéT pacioJ1IOKeHnn 
WM WonepedHplx ceyeHnit NOABOXHBIX KOHCTpyKUHH B HOTOKe KpHTH4eCKOLO PeKHMAa 
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AHHOTalna 

Beedenue. J\nmuusie u y3Kve B MOMepeyHHKe KOHCTpyKUMM MOPCKHX IHEprosOOKIBalOWHX CHCTeEM HaxXOATCA MO T0- 
CTOAHHBIM BO3eMCTBHeM TeyeHH wu BOTH. [uapodqMHaMMyecKne Harpy3KH ABJIAIOTCA Pe3yJIbTATOM B3aMMOJeHCTBUA 
MOABOAHBIX TpyOOMpoBOOB, WaHToKaOeseH, oop OOOpyAOBaHHA C MOTOKOM 2KMAKOCTH HM MPHBOAAT K OOpa3z0BaHHto 
BUXpel B 30He 3a KOHCTpyKuMaAMu. BuxpeoOpa30BaTeJIbHble CHJIbI CILYKAT HCTOUHHKOM UMKIMYeCKOrO Harpy2KeHHA UH 
TIOCTEMCHHO YCKOpAIOT YCTaIOCTHOe pa3spylIeHHe, YTO MOXeT MpHBeCTH K aBapHAM. OJHUM U3 CHOCOOOB CHH>KeHHA 
Harpy30K Ha NOABOJHbIC KOHCTPYKUMH ABJIACTCA H3MCHEHHe (POPMBI HX MOMepe4yHoro CeyeHHA C YUCTOM pexKHMa HOTOKAa. 
HeyoctaTouHo u3y4eHo, KaKHM 06pa30M HTOTOBbIe rH pOMHAaMMYeCKHe Harpy3KU 3aBHCAT OT (POPMBI WoMepeyHoro ce- 
YeHHA HW B3AMMHOTO paciOJIOXKeCHHA HA3BAHHBIX BbILIe 3ICEMCHTOB CHCTeM, HaXOJ{ALIMXCA B paBHOMEPHOM KPHTH4YeCKOM 
mioToke. IIpeyctapyeHHad Hay4ynad padoTa MpH3BaHa BOCHOJHUTS ITOT Mpoben. Lemb ucceqyoBaHHA — PpaccMOTPeTb B 
JJaHHOM KOHTeKCTe 3Ha4eHHe paccTOAHHA MeX*KAy KOHCTpyKUMAMH, a TakKxKe HaMuve NosyKpyrion D-oOpa3sHoi KOH- 
CTPpyKIUMH, pa3sMeléHHOH Nepey rpynioli 43 Tpéex WWM APOB C pasHbIMM MonepeuHbIMM CeyeHHAMH. 

Mamepuanoi u memooot. Jna wicneHHoro MoyeuMpoBaHHA BUXpeoOpa3z0BaTeJIbHbIX CHI HCIHOIb3OBaJICA MCTO], BbIMHC- 
JMTeIbHOM JHHAaMUKH (pIOMAOB B IporpamMMe ANSYS Fluent aa uA ApoB AMametpom D = 0,3 m. MoyemupoBanue 
BBIMOJIHCHO Me€TOJOM HelIpHcoeqHHéHHBIX BUxpeli DES, koTopblii coueTaeT B cebe NpeMMylecTBa MeTOa ycpeqHeH- 
Horo 10 PeitHombacy ypaBHenua Happe-Ctoxca RANS nu metoyia kpynHbIx Buxpeli LES. B kayectBe oOnexta uccieqo- 
BaHHA pacCMaTpHBallach CHCTeMa, COCTOAMAd W3 YeTLIPEX KOHCTPYKIMM B BbIYMCIMTeIbHOM JOMeHe B 2D, BKuIOUAA 
CTOALIMM BbIMe 10 TEYeEHHIO NOMYKpyribli WWIMHAp HU OCHOBHY!O rpyninly 43 Tpéx WHIMHApOB Kpyrioi, KBaypaTHoH u 
POMOoOBH HOM (POpMbI WoMepeyHoro CceyeHHA. ITH KOHCTPYKUMU B yCIOBHAX H€YCTAaHOBHBIIeroca polecca HaxOWATCA 
lO AelicTBHeM paBHOMepHOro NoTOKa KpuTMYecKOro pexkuma pH Re = 2,5x10°. 

Pe3yibmamot ucciedosanua. B pe3yibTaTe MOJeIMpOBaHHA MOJYYCHI MATb HAOOPOB AaHHBbIX JIA M3MCHAIOLIMXCA BO 
BpeMeHH KOOPPUUMEHTOB BUXpeoOpa3z0BaTeNbHBIX MOAbEMHOM CHJIbI HM CHJIbI COMPOTHBIICHHA: JIA OCHOBHOM CHCTeMBI 
M3 MOYKpyrsiol, Kpyriol, KBadpaTHOM H pOMOOBHAHOM KOHCTpyKUMH, a TakoKe JIA YeTbIPEX CHCTEM 3 TOJIbKO TIOJy- 
KPYTJIBIX, TOJIBKO KPYTJIBIX, TOJIBKO KBaPaTHBIX H TOJIbKO pOMOOBHHBIX KOHCTpyKUHH. JJOMOJHUTeIbHO IpOBeyGH aHa- 
JIM3 BIMAHHA paCCTOAHHA MOK Ly KOHCTPyYKUMAMH Ha aMIWVIMTyAy KoeOaHnit KoIpPUUHeHTOB rH pOAMHAaMHYeCKHKX CHIL. 
Tlomy4enHble pe3yIbTaTbl MpeACTaBIeHbI B BUC KOIPPUUMEHTOB NOAEMHOM CHJIbI M CHJIbI COMpPOTHBIICHHA B MHa- 
MHKe, aHasI43a YaCTOT HM KOHTYPOB MOJIeH CKOPOCTH, WaBIeHHA, 3aBUXpEHHOCTH. Pe3ybTAaTbI NO3BOJIAKOT YCTAHOBUTb 
HOJIOXKUTEIbHOE BIMAHME CTOALMICH BbIMEe MO TEYCHHIO NOUYKpyrsiol KOHCTpyKWHH Ha CHW KeHHE CHJIbI CONPOTHBIICHHA 
Ha WeHTpasIbHy1O KOHCTPyKIHH0 B rpylitie 43 TpéxX WHIMHApOB HDKe 110 TeYeHHIO. 

O6cystcoenue u 3akitovenue, Pe3yibTaTbI IpOBeCAEHHBIX HUCCIeTOBaHHM MO3BONAIOT IPHHAMaTb OOOCHOBAHHEIEe pellie- 
HUA [IX PaCCTAaHOBKM MOPCKHX KOHCTpyKUMM B pyle 43 YeTLIPEX OOLEKTOB B 3ABHCHMOCTH OT (POpMBI HoMepeyHoro 
ce4ueHHaA HU paccTOAHMA MexKy HMMM. YcTaHOBKa NOYKpyrioH KOHCTPyKUMH BIE TO TEYeHHIO NO3BOJIAeT CHH3HTB 
TH pOAMHAMHYeCKy!O CHJIY COMPOTHBJICHHA Ha WeHTpasIbHy!O KOHCTPyKUMHO B rpyniie 43 TpéxX KOHCTpYKUMH HWKe 110 
TeC4eHHIO, YTO 3aMeIAeT C€ YCTANOCTHOE paspylleHHe HW YBeJIMYMBAeT CPOK SKCIIyaTauHu. 
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Introduction. Operation and construction of modern offshore systems, specializing on the energy production, 
extraction of resources or carbon capture and storage, require evaluation of the impact of environmental flows on 
equipment and structures. An increased fatigue failure in subsea structures, such as pipelines, risers, cables, piles, 
equipment supports, may come from the vortex shedding phenomenon. The problem is particularly important when 
slender structures are designed to reach deep waters to connect subsystems together. The layout of subsea systems comes 
with the arrangement of structures with different geometry, hydrodynamic properties, and their position in proximity to 
each other. The interference of wakes from these structures and vortex formation patterns is sometimes challenging to 
predict at very high Reynolds numbers due to the turbulent nature of the flow. 

Differences in the flow over a standalone cylinder and two cylinders in tandem are discussed in [1], and three vortex 
shedding regimes are identified for tandem structures. These vortex formation modes include the extended-body regime 
at 1.0 < L/D < 1.8, where L/D is spacing ratio, commonly used to quantify the distance between centres of neighbour 
structures. So that, Z corresponds to the distance, and D is the diameter of the structure. In [1], increase of the spacing 
ratio to 1.8 < L/D < 3.8 leads to the reattachment regime, where shear layers detach from the upstream structure and 
reattach to the front side of the downstream structure, so that vortices are formed behind this downstream object. Further 
growth of the spacing ratio, above L/D > 3.8, introduces the co-shedding regime, where a separate vortex is formed from 
the upstream structure and from the downstream structure. Another fundamental research investigated vortex dynamics 
in details through experimental research [2]. One of the following fundamental studies [3] experimentally investigated 
the vortex shedding frequencies of two staggered identical circular cylinders with the Reynolds number Re varying from 
3.2x10* to 7.4x104, and two fixed side-by-side cylinders at the Reynolds number of 2.5x10* were earlier considered in [4]. 
These investigations provide an important foundation for modern studies in terms of the known effects in fluid forcing 
and vortex shedding patterns. 

Significantly more recent investigations are performed numerically [5] to study the effect of spacing on loads and 
vibrations for two tandem cylinders at subcritical Reynolds numbers, and for specific cases, like a group of mixed large 
and smaller structures [6]. The latter work [6] numerically investigates fluid force coefficients and observes the vortex 
formation pattern on three identical rigid circular cylinders in proximity to a square cylinder. A parametric study is 
conducted in [7] for three identical stationary circular and D-shaped cylinders placed close to a square cylinder at 
Reynolds number 3900 in both linearly and parabolic sheared flows. 

Considering the impact of cross-sectional shapes further, a numerical study is conducted by [8] for a flow over six 
identical stationary cylinders having different cross-sectional shapes at Reynolds number of 2.5x10° in the uniform and 
linearly sheared flow. Rectangular cylinders are investigated in details in [9, 10, 11], where one of the most impactful 
factors for hydrodynamic loads is the aspect ratio of rectangle sides. The works [9, 10] provide new experimental data 
and attempt to develop semi-empirical methods of predicting the response of structures. Further steps in improving the 
modelling approaches for the structural vibration of rectangle-shaped objects under the hydrodynamic excitation are 
performed in [11]. Another branch of studies considers a flow over a sphere [12, 13, 14], while still leaving issues of the 
impact of the cross-sectional shape of subsea structures open. 

Further research on diverse cross-sectional shapes is performed in [15, 16, 17], where triangular and diamond cross- 
sections are studied in comparison. Research [15] is focused on the sensitivity to the corner sharpness for the diamond 
(rhomb) cross-section. Work [16] investigates effects from diverse cross-sections for the system with a rotational degree 
of freedom, when subjected to flow-induced vibration, and study [17] uses cross-sectional shapes for the energy 
harvesting with fluid-structure interaction. 

Following published results, the current work aims to investigate the drag reduction when three structures of different 
cross-sectional shapes are located around a circular cylinder to observe the wake interference and the vortex formation 
pattern between these structures for the spacing ratio L/D varying from 1.67 to 2.83 in the uniform flow at the Reynolds 
number of 2.5x10° with the computational fluid dynamics approach. Specifically, the drag reducing ability of the upstream 
structure is of the research interest. The considered layout is a combination of cylinders placed in tandem, side-by-side 
and staggered position, with a mixture of cross-sectional shapes. 

Materials and Methods. A system of four cylinders with a diameter (characteristic size) of 0.3 m with different cross- 
sectional types is considered in this study. The selected cross-sectional shapes and positions of structures are shown in 
Figure 1. Fluid forces and the flow interference are studied for the spacing ratio of L/D varying from 1.67 to 2.83 in the 
uniform current at the Reynolds number of 2.5x10° and the corresponding velocity of 0.837 m/s. The study focuses on 
the specific layout, where structures are positioned relative to each other in a mixed tandem (cylinders | and 3), side-by- 
side (cylinders 2 and 3, and cylinders 3 and 4) and staggered (cylinders | and 2, and cylinders | and 4) configuration. 
Cross-sectional shapes considered include half-circle, square, circle, and diamond. 
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Periodic 


Shadow 


Fig. 1. Considered structures in the computational fluid domain 


Computational fluid dynamic simulations are performed in the domain with a size of 47D x 27D. The top and bottom 
boundaries are located at a distance of 13.5D away from the center of the circular cylinder, periodic and shadow properties 
are assigned to these boundaries. The left boundary serves as a velocity-inlet, which is located at the distance of 20D from 
the centre of the circular structure 3 in the domain. The value of gauge pressure is set to zero at the pressure-outlet set at 
the right boundary. No-slip conditions are applied to the cylinders. 

The flow around cylinders is simulated using the computational fluid dynamics software ANSYS Fluent, where the 
finite volume method is implemented to solve the Navier-Stokes system. The incompressible flow is considered, and the 
2D DES transient simulations are conducted with the A-w SST turbulence model. Time integration is performed using 
the second-order implicit transient formulation with a time step of 0.01 s, and the PISO algorithm is used as the solver. 

The DES approach connects capabilities of the Reynolds-Averaged Navier-Stokes (RANS) and Large Eddy 
Simulation (LES) methods [8]. The RANS governing equations for the incompressible flow are as follows: 


8(pu; 
= =0, (1) 
A(pu;) 9 ,— —— _ op , Oy p 
A #3 (Puie tenia = <P Ox) (2) 


where p is mean pressure, u, is average Cartesian components of the velocity vector, pu;u; are Reynolds stresses, p is 


density of the fluid, and t,; is mean viscous stress vector components, which could be expressed as: 


= OU 
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where pl is dynamic viscosity. 
The Large Eddy Simulation (LES) system of equations for the incompressible flow can be written in the following way: 
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where U; and P represent the resolved filtered velocity and pressure, respectively. 
The diffusion term of the DES model is given by 
Y, =pB'koF pes, (6) 
where " is a constant, k stands for fluctuation of the turbulent kinetic energy, w is specific energy dissipation rate, and 
F pgs is as follows: 


Foes = mas a, ‘} (7) 


des max 


Advanced Engineering Research (Rostov-on-Don). 2024;24(2):135-147. eISSN 2687-1653 


where Cas is a constant, Ajax is local maximum grid map A = (Aj, Ao, A3)!*. Further, L; is turbulent length scale: 


Je 


= oe 8 
(Bo (8) 
The DES-SST model uses the following zonal formulation: 
L 
Foes = mas| : (1 Fsr).t} (9) 
des‘ max 


where F'ssr = 0, F1, Fo, and F), F are mixed functions of the SST model. 

Table 1 below provides results of the mesh independence test for the uniform flow of Re = 2.5105. The mesh settings 
are adopted from [7], and the accuracy of the grid is demonstrated by comparisons in Table 1. All subsequent analysis in 
this paper is performed with Mesh 2, and the results include signals and frequencies of the fluid force coefficients and an 
indication of the vortex shedding pattern. The drag force fluctuations are presented in this work in terms of the drag force 
coefficient Cp which comprises the mean drag coefficient Cpo and the fluctuating drag coefficient Cp", as follows: 


Che Cages. (10) 


The lift force fluctuations are presented using the lift coefficient C;. 


Table 1 
Mesh independence test results 
Cases Cbo Number of cells Strouhal number 
Re = 2.5x10° 
Current study 
Mesh 1 0.98 63,345 0.24 
Mesh 2 1.08 86,478 0.24 
Mesh 3 1.08 131,041 0.24 
Published data 
Lehmkuhl, et al. (2014) (LES) [18] 0.833 - 0.238 
Achenbach&Heinecke (1981) (Experiment) [19] 1.135 - 0.230 
Re = 3,900 
Current study, Mesh 2 0.93 86,478 0.18 
Wornom, et al. (2011) (VMS-LES) [20] 0.99 = - 
Re = 3.6x10° 
Current study, Mesh 2 0.4100 86,478 = 
Porteous, et al. (2015) (URANS) [21] 0.4206 - =: 
Nazvanova, et al. (2022) (URANS) [22] 0.4657 74,496 = 


Results. In this study, simulations are performed in two series. The first series is focused on recognising the overall 
effect of various cross-sectional shapes, placed at L/D = 2.00 from each other. Cylinder numbers here correspond to the 
ones used in Figure 1, according to the cylinders position. The calculation results for this set are reported in Table 2, in 
comparison to the case of structures with mixed cross-sections. The second series provides an insight into the impact of 
L/D ratio on hydrodynamic loads observed for the mixed cross-section case only, as in Figure 1. These results are 
summarised in Table 3 and allow defining the drag reducing effect of the D-shape upstream structure on loads when 
cross-sections are different. 
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Simulation results for the same arrangement with different cross-sectional shapes for L/D = 2.00 


Table 2 


LID =2.0 ei oi nie ia All circular | All square All D-shaped All diamond-shaped 
Groseeectons in Riga structures structures structures structures 
Cylinder 1 
Coo 0.45 0.28 0.49 0.48 0.46 
ch 0.19 0.14 0.43 0.28 0.17 
CL 0.18 0.04 0.07 0.13 0.28 
Cylinder 2 
Coo 0.87 0.39 0.89 0.64 0.70 
cf 1.03 0.38 1.38 0.49 0.60 
Cr 1.38 0.53 0.07 0.11 0.07 
Cylinder 3 
Co 0.3 0.24 0.26 0.35 0.70 
cE 0.45 0.37 0.91 0.68 0.63 
CL 1.03 0.73 1.26 0.30 0.73 
Cylinder 4 
Cpo 0.98 0.39 0.83 0.61 0.99 
cf 0.90 0.27 0.96 0.44 1.05 
Cr 1.22 0.81 1.63 0.37 1.11 
Table 3 
Simulation results for the mixed cross-sections at various L/D ratio 
Cylinder 1 Cylinder 2 Cylinder 3 Cylinder 4 
L/D 
Co | CH CL Co | Ch CL Co | Ch Cr Co | Ch Cr 
1.67 0.14 0.18 0.15 0.98 1.17 1.44 0.48 0.57 0.66 0.93 0.75 1.13 
1.83 0.33 0.21 0.18 0.93 1.54 1.40 0.33 0.50 0.86 1.01 1.36 1.15 
2.00 0.45 0.19 0.18 0.87 1.03 1.38 0.30 0.45 1.03 0.98 0.90 1.22 
2.17 0.49 0.05 0.19 0.81 0.14 1.11 0.25 0.06 0.83 0.96 0.38 1.10 
2.33 0.37 0.12 0.00 0.78 0.34 1.03 0.25 0.06 0.28 1.00 0.07 0.79 
2.50 0.57 0.16 0.20 0.76 0.89 1.18 0.26 0.47 0.77 1.01 0.64 1.06 
2.67 0.57 0.02 0.27 0.78 0.16 1.37 0.28 0.08 0.82 1.02 0.03 0.96 
2.83 0.58 0.21 0.21 0.84 0.91 1.70 0.20 0.52 1.16 0.93 0.68 1.03 


Comparison of all circular, all square, mixed (as in Fig. 1), all diamond and all D-shapes with each other in Table 2 
reveals relatively lower hydrodynamic loads for all four structures observed for the circular shapes at the L/D = 2.00. The 
largest mean drag coefficient here is experienced by the fourth structure in both mixed-shaped and all diamond-shaped 
arrangement. The highest maximum fluctuating drag coefficient of 1.38 is observed for cylinder 2 in the all square-shaped 
arrangement. The largest maximum amplitude of the lift coefficient of 1.63 also belongs to the all square-shaped 
arrangement, but corresponds to cylinder 4. Mixed-shaped arrangement (or basic case of structures with alternate cross- 
sections, as in Fig. 1) at L/D = 2.00 demonstrates a relative consistency in large amplitudes of the lift coefficient for 
cylinders 2, 3, 4, which makes estimations of hydrodynamic loads for this arrangement more important, due to higher 
expected loads. 
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Fig. 2. Time histories of fluctuating drag coefficients for four cylinders of a different cross-sectional shape with L/D = 2.0: 
a — cylinder 1; b — cylinder 2; c — cylinder 3; d— cylinder 4 
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Fig. 3. Time histories of lift coefficients for four cylinders of a different cross-sectional shape with L/D = 2.0: 
a — cylinder 1; b — cylinder 2; c — cylinder 3; d — cylinder 4 


Further observation of signals of the fluctuating drag and lift coefficient in Figures 2 and 3 reveals meaningful instabilities 
in forces experienced by all structures with square shapes, and for structures 2, 3, 4 with diamond and mixed shapes. The 
comparison shows that circular and D-shaped structures would experience lower and more stable fluid loads in this arrangement. 
Table 2 and Figures 2 and 3 confirm that cylinder 1 in the shape of a circle or D-shape provides reducing fluid loads for the 
three downstream cylinders. This substantiates the common interest in further exploration of the effect of the D-shaped structure 
on reducing the drag force in the arrangement with all alternative cross-sectional shapes. 
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This effect is studied in details in the next (second) simulation series, presented for each cylinder in Figures 4—7 in 
terms of the fluid forces and in Figures 8, 9 — in terms of the fluid flow characteristics for the considered computational 
domain. Figure 4 illustrates fluid loads on the upstream structure, where the most unstable signal (at L/D = 2.17) belongs 
to fluctuations of the drag force. Figure 4 c also indicates presence of multiple frequencies in signals of the fluctuating 
drag coefficient, while a single dominating frequency can be identified for the lift force coefficient in Figure 4 d. 
Figures 4 c and 4 d demonstrate that the frequency of both lift and drag forces generally increases with the growing L/D 
ratio, and the maximum frequency is indicated by signals at L/D = 2.67 and 2.83. 
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Fig. 4. Time histories of fluid force coefficients for cylinder 1 in the uniform flow: 
a — fluctuating drag coefficient; b — lift coefficient; c — fluctuating drag coefficient FFT; d — lift coefficient FFT 


Relatively similar complexity of frequencies of the lift force is observed in Figure 6 d for cylinder 3, where one to 
two dominant frequencies could be clearly identified. At the same time, more than two frequencies are observed in 
Figure 6 c for each signal of the fluctuating drag force. The pattern of growth in the overall dominant frequency with the 
increasing L/D is still recognisable for cylinder 3, similar to cylinder 1. Combination of frequencies is even more complex 
for cylinders 2 and 4, as comes from Figures 5 c—d and 7 c-d, this is not possible to indicate clear dependences in the 
frequencies of the fluctuating drag force. Some resemblance of the found growth trend could be still observed in Figure 5 d 
for the lift force coefficient for cylinder 2. This provides evidence for the generally unstable nature of hydrodynamic 
loads acting on the square cross-section shown in Figure 5. 

Impact of the L/D ratio on the mean drag coefficient for cylinder 1, based on Table 3, is partially confirmed: apart 
from a couple of deviations, the mean drag force increases with the growth of L/D. Data for cylinders 2 and 4 do not 
indicate a specific pattern, as the values fluctuate back and forth within 10% from the initial mean drag coefficient at the 
smallest L/D. Cylinder 3, on the contrary, indicates a stable pattern of the reduced mean drag coefficient with the increased 
L/D, so that the reducing ability of the upstream D-shaped structure is evident, but for the central cylinder only. The 
largest mean drag coefficient of 1.02 is linked to cylinder 4 at L/D = 2.67. The considered range of L/D allows observing 
an important transition from the strong to minor interference in the wake of the three paired structures. 

The highest fluctuating drag coefficients of 1.54 and 1.36 are linked to cylinders 2 and 4, respectively, both observed 
at L/D = 1.83. The feature of the maximum amplitude of the fluctuating drag coefficient, indicated in Table 2, is in absence 
of a specific dependence from the L/D, the values rapidly change from near zero to relatively high with a small increment 
of change in the ratio. The fluctuating drag coefficient has generally the lowest amplitudes for cylinder 1, average 
amplitudes — for cylinder 3, and the largest amplitudes — for cylinders 2 and 4. 

The largest maximum amplitude of the lift coefficient occurs at L/D = 2.83 for cylinder 2. The lift coefficient 
generally resembles the distribution, similar to the maximum amplitude found for the fluctuating drag coefficient: the lift 
force appears to be the smallest for cylinder 1, relatively intermediate — for cylinder 3, and the highest — for cylinders 
2 and 4, with no specific pattern linked to the Z/D increase and reduction. This allows us to conclude that the ability to 
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reduce hydrodynamic forces by placing the upstream D-shaped structure in front of the array is limited. The force 
reduction is observed mainly for the structure placed in tandem downstream, and the effect is most pronounced for the 
mean drag coefficient, with some reduced effects also seen for the fluctuating drag and lift coefficients. 
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Fig. 5. Time histories of fluid force coefficients for cylinder 2 in the uniform flow: 


fl 
D 


Fluctuating drag coefficient, 


106 107 108 109 110 
Flow time, t (s) 


a) 


So 
we 


Power spectrum density, PSD 


a — fluctuating drag coefficient; b — lift coefficient; c — fluctuating drag coefficient FFT; d — lift coefficient FF 


0.4 0.6 0.8 1.0 1.2 1.4 
Frequency, f (Hz) 
°) 


Lift coefficient, C, 


0.8 : 
105 106 107 108 109 110 
Flow time, t (s) 
b) 

14 T T 
—— L/D=1.67 
12} — LD = 1.83 
L/D = 2.00 
L L/D =2.17 

10 

—— L/D =2.33 


Power spectrum density, PSD 


0.50 0.55 0.60 0.65 0.70 0.75 0.80 


8 ——— LID = 2.50 


0.3 0.4 0.5 0.6 0.7 


L/D = 2.67 
L/D = 2.83 


Frequency, f (Hz) 
d) 


Fig. 6. Time histories of fluid force coefficients for cylinder 3 in the uniform flow: 
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Fig. 7. Time histories of fluid force coefficients for cylinder 4 in the uniform flow: 
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Fig. 8. Contours of the flow characteristics for L/D = 1.67 at 200 seconds: 
a — velocity contour; b — vorticity contour; c — pressure contour; d — streamline 
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Fig. 9. Contours of the flow characteristics for L/D = 2.83 at 200 seconds: 
a — velocity contour; b — vorticity contour; c — pressure contour; d — streamline 


Figures 8—9 show the velocity, vorticity, pressure and streamlines of the flow around cylinders for some selected L/D, 
where both proximity and wake interference among the cylinders are presented for the time step of 200 s. The flow around 
cylinders is complex, and vortex formation patterns are highly affected, as the distance between the cylinders increases. 
The proximity interference is observed for cylinders 2, 3, and 4, alternate single vortices are shed on the downstream side 
of these structures. For the wake interference at L/D = 1.67, free shear layers separate from the upstream cylinder | and 
reattach themselves to the upstream side of cylinder 3, and a vortex street is only formed at the downstream side of 
cylinder 3. At this distance, a broad region of wake is created at the downstream side of cylinders 2, 3, and 4. As the L/D 
increases above 2.00, there are vortices formed at the upstream side of cylinder 3, in the wake of cylinder 1. A vortex 
street is also formed at the downstream side of cylinder 3, with formation of 2S vortices. Figures 8 a, 9 c demonstrate a 
group of minor vortices formed following the diamond-shaped cylinder 4 and a vortex pair formed in a similar to S+P 
vortex cycle past cylinder 2. 

Discussion and Conclusion. The 2D numerical simulations are performed in this work for cylinders with different 
cross-sectional shapes at the Reynolds number of 2.5x10° using the DES approach. The considered cylinders are studied 
in a complex position of an upstream D-shaped structure in front of three paired structures, with the aim to investigate the 
drag-reducing ability of this specific layout, observe the flow complexity, the wake interference from each structure, and 
vortex formation patterns. 

The following conclusions could be made from this study: 

1. The ability to reduce hydrodynamic forces by placing the upstream D-shaped structure is mainly limited to the 
structure placed in tandem downstream, and the effect is most pronounced for the mean drag coefficient. 

2. Overall, the mean drag coefficient of cylinders is observed to be affected by varying L/D, with the main effect on the 
mean drag coefficient of cylinder 1, which grows with increasing L/D, and of cylinder 3, which reduces with increasing L/D. 

3. Competition of frequencies is observed for the fluctuating drag coefficient for all structures and lift coefficient 
signal of cylinders 2 and 4. This competition is due to the joint effects of both the uniform current and wake interference, 
which intensifies at a lower L/D in terms of changes to the resulting vortex street. 

4. Both proximity and wake interference among the cylinders are observed. The flow around cylinders is complex, 
and vortex formation patterns are highly affected as the distance between the cylinders increased with 2S being the major 
vortex type formed and shedding the additional vortices from the square and diamond structures. 
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The study generally contributes to the field of knowledge by advancing our understanding of fluid-structure 
interactions, drag reduction strategies, and vortex dynamics, with potential applications in offshore energy systems. The 
current work contributes to the development of the drag reduction strategies through analyzing the impact of the upstream 
D-shaped structure on downstream cylinders. Understanding how different structural configurations affect drag can 
inform the design of more efficient systems in various engineering applications, such as subsea transportation of fluids. 
By observing flow complexity, wake interference, and vortex formation patterns, this study contributes to the 
understanding of fluid dynamics around complex geometries. This knowledge is crucial for optimizing the performance 
of structures in environments where fluid flow plays a significant role, such as subsea engineering. The results of the 
present research highlight the effect of varying the aspect ratio L/D on drag coefficients to inform engineering designs of 
similar arrangements. This study reveals the intricate vortex dynamics and shedding patterns, particularly concerning the 
proximity and wake interference among cylinders. Understanding these phenomena can aid in predicting and controlling 
flow behavior around complex configurations, leading to more efficient designs and better performance in practical 
applications in offshore systems. 
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Abstract 

Introduction. Devices for collecting and storing energy from the external environment are low-power sources of electric 
energy that are actively used. The autonomous devices for monitoring the damaged condition of various structures include 
them as well. The working element of these devices is a piezoelectric generator (PEG) — a converter of mechanical 
energy into electrical energy. The design of PEG is associated with the preliminary construction of their mathematical 
and computer models, with the help of which the calculation and optimization of structures is carried out. One of the ways 
to model and calculate PEG is to develop approximate calculation methods based on applied theories. The applied theories 
for calculating bending vibrations of multilayer piezoactive plates are known and previously developed in the literature. 
However, in the scientific literature there is not enough information about bending and shear vibrations as a tool for 
improving the efficiency of engineering calculations of the described structures. The objective of this work was to develop 
an applied method for calculating bending and shear vibrations of piezoceramic plates, including porous ones. 

Materials and Methods. Piezoceramics PZT-4, including porous ones, were used as the piezoactive material of the plate. 
When using porous ceramics, the rigidity of the structure decreased to a greater extent than the piezoelectric modules, 
which made it possible to obtain a more effective PEG under mechanical action. The mathematical formulation was 
carried out within the framework of the linear theory of electroelasticity with plate polarization in thickness. The sides of 
the plate were electrodated, the right side was fixed, and a smooth contact in the vertical wall was set on the left side. 
Steady-state vibrations of the plate were caused by pressure on the front surfaces of the plate or the difference in electrical 
potentials at the electrodes. To calculate the characteristics of PEG, the authors proposed an applied theory based on 
hypotheses about the distribution of characteristics of the stress-strain state and the electric field. 

Results. Transverse vibrations of a piezoceramic plate in the low-frequency region (below the first bending-shear 
resonance) were studied. Due to the fact that the mathematical formulation was considered within the framework of the 
linear theory of elasticity, the problem was divided into the sum of two. The first one took into account the mechanical 
effect: a distributed load and a transverse force at the left end acted on the front surfaces of the plate, and the potentials 
at the electrodes were zero. In the second task, there were no mechanical loads, but the potential difference was set at the 
electrodes. Based on hypotheses about the distribution of deformations, mechanical stresses and electric potential, both 
problems were reduced to a system of ordinary differential equations and boundary conditions. Comparison with the 
results of calculations by the finite element method in the ACELAN package showed the adequacy of the proposed applied 
theory in the low-frequency region. 

Discussion and Conclusion. Since the formulation of the problem was considered in the linear theory of electroelasticity, 
and the low-frequency region was studied, the work succeeded in dividing the problem of bending-shear vibrations of a 
porous piezoceramic plate into two: bending — with mechanical action at zero potentials, and shear — when setting the 
potential difference and zero mechanical action. The corresponding hypotheses about bending and shear were used. Two 
systems of ordinary differential equations and boundary conditions, which were solved analytically without the use 
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of “heavy” finite element packages, were constructed. To compare the results and confirm the adequacy of the proposed 
method, the finite element modeling of such tasks was carried out in a specialized ACELAN package. The comparison 
showed that the error in determining displacements and electric potential when using this approach, in the case of setting 
mechanical loads and potential differences, did not exceed 6%. The method developed in the paper can be applied in the 
design of piezoelectric generators for energy storage in the low-frequency region. 


Keywords: energy collection device, piezoelectric generator, porous ceramics, plate bending, plate shear, applied theory 
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AHHOTalna 

Beedenue. YctpolictBa cOopa u HakOMJIeHHA JHEP M3 BHEeIMIHeM CpebI MpeCTaBIAIOT COOOM MasIOMOMIHbIe HCTOU- 
HHUKH 3JICKTPHYeCKOM IHEPrHu, KOTOPbIe AKTHBHO HCIOMb3YIOTCA, B TOM YHCIIe B ABTOHOMHBIX IpHOopax MOHHTOpHHra 
MOBPeK CHHOTO COCTOAHHA pasIM4YHbIX KOHCTpyKUMH. PaOouuM 3IeEMeHTOM 3THX YCTPOHCTB ABIIACTCA IIbe309IEKTpH4e- 
ckui renepatop (IISI’) — mpeoOpa3oBatemb MexaHH4eckol 9HeprHuu B WIeKTpuyeckyro. Konctpyuposaune IIDT caa- 
3aHO C IIpeABapHTesIbHbIM MOCTpOeCHHeM HX MaTeMaTHYeCKHX HM KOMIBIOTEPHbIX MOeJeH, C MOMOLI[bIO KOTOPbIX IPOH3- 
BOJHTCA pacueT HU ONTHMU3AIIMA KOHCTpyKuMi. OHUM v3 crocoboB MoyemupoBaHua u pacueta IIOT saBnaetca pa3pa- 
OoTKa MIpHOMWKeHHBIX METOJOB pacueTa Ha OCHOBe IPUKIaTHbIx Teopuii. B mHTepaType U3BeCTHBI HM paHee pa3spa0oTaHbl 
IIpHKaHble TeEOpH pacueTa H3rHOHBIX KOIeOAaHHM MHOFOCOMHBIX Mbe3OaKTHBHBIX IWiacTHH. OqHako HHpopMalnH 06 
W3PHOHO-CABHTOBBIX KOeOaHHAX, KAK HHCTPyMeHTe MOBbILIeHHA I:pPCKTHBHOCTH HHKCHEPHBIX PaCUeTOB OMMCaHHBIX 
KOHCTpyKUMH, B HaydHOH MTeparype HegocTaTouHo. Llembio HacTosUel paOoTs! ABIIAIach paspaboTKa IpHkKaqHoro 
MeTOJja pacueTa H3rHOHBIX MH CABHTOBBIX KoJeOaHHi be3OKepaMHUeCKHX IIaCTHH, B TOM YMCIIe MOPHCTEIX. 
Mamepuanoi u memoooi. B xauectTBe 1be30aKTMBHOrO MaTepHalia MWaCTHHbI MCHOb3yeTCA Wbe3s0Kepamuka PZT-4, B 
TOM 4ucIe Mopuctas. [pu ucnoub30BaHHH MOpHCTOM KepaMUKH %KCCTKOCTb KOHCTpYKUMH YMeHbIUWaeTCA B OObIUeH CTe- 
TICHH, UCM MIbe3OMOJYIIN, YTO MO3BOIACT MOMYYNTh Ooee sspexTHBHBIM ITST mpu Mexanwyeckom BosyelicTBuu. Mate- 
MaTH4eCKad MOCTaHOBKa OCYIIeCTBIeCHa B PaMKaX JIMHeCHHOM TeEOPHU ZICKTPOYNPYFOcTH Ip NOMAPH3alMH TWIAaCTHHbI 
110 ToMMHHe. boKOBbIe CTOPOHBI MJIaCTHHI IIEKTPOAMPOBaHbl, IpaBad CTOPOHa 3akpelieHa, a Ha JIeBOn 3ayaH riayKn 
KOHTAKT B BEPTHKaJIbHOM CTeHKe. YCTaHOBMBIIHeCA KOJICOaAHHA IIaCTHHBI BbI3IBAIOTCA aBJICHHeM Ha JIMI[eBbIe IOBEpX- 
HOCTH IJIlaCTHHbI WIM pa3sHOCTbIO 3ICKTPHYeCKHX MOTCHUMAIOB Ha WIeKTponax. [1a pacueta xapaxtepucTuK IID B pa- 
OoTe MpeaaraeTca MpHKayqHad TeopuA, OCHOBaHHad Ha PHMOTesax O pacnipeseeHHU XapakTepHCTHK Halips»KeHHO-]e- 
(POpMUpOBaHHOTO COCTOAHHA HM 3IEKTPHYECKOTO MOA. 

Pe3yivmamot ucciedoeanua. PaccMoTpeHbl MonmepedHble KoIeOaHHA Ibe30KepaMH4eCKON MJIaCTHHbI B HA3KOUaCTOTHOM 
oOsacTH (HYKe TepBoro H3rMOHO-CABHTOBOrO pe3zoHaHca). B cusy Toro, YTO MaTeMaTHUeCKad MOCTaHOBKa paccMOTpeHa 
B paMKax JIMHeMHOM Teopuu yimpyrocTnH, 3ayjaya pa3zyeuacb Ha cyMMy Byx. B nepBolt yaHTEIBaIOCb MexaHH4eckoe 
BO3elCTBHe: Ha JIMI[eBble MOBEPXHOCTH IIaCTHHbI JelcTByeT pactipeyeseHHad Harpy3ka MU MomepedHas cia Ha JIEBOM 
KOHIe, a MOTCHUMasIbI Ha IICKTposax paBHbl HyO. Bo BTopo 3ayave MexaHH4eckHe Harpy3KH OTCYTCTBOBAJIH, HO 3a- 
JjaBalacb pa3HOCTb MOTeHUMaOB Ha 9eKTpopax. Ha ocHoBe rumoTe3 0 paciipeyzereHuu AepopMalHi, MexaHHyeckux 
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HallpsoKeHHH UM 9IEKTpH4eCKOrO MOTeHMasa OOe 3aa4H ObIM CBECHBI K CHCTeMe OOBIKHOBeHHBIX JHPepeHuMasIbHBIX 
ypaBHeHuli U TpaHW4HBIX yCOBMH. CpaBHeHue C pe3yJIbTaTaMM pacueTOB MeTOJOM KOHC4UHBIX 3JIEMCHTOB B TlakeTe 
ACELAN noka3asIv ayjeKBaTHOCTb MpesOKeCHHOM NpHKIaqHON TeOpHu B HU3KOUAaCTOTHOM OOacTH. 

O6cystcoenue u 3aknroyuenue. TlockolbkKy NOCTaHOBKa 3afa4 pacCMaTpHBasaCcb B JIMHEMHOM TeEOpHU BICKTPOyMpyrocTu U 
W3y4alacb HH3KOYAaCTOTHAaA OOACTb, B padoTe yAaslocb 3ayauy OO H3THOHBIX HM CABHTOBbIX KOJIeOaHHAX IWIaCTHHBI 13 T10- 
PUHCTOH Mbe30KepaMUKH pa3z{eIMTb Ha Be: H3rMOHy!O — C MeXaHHYeCKHM BO3eHCTBHeM IIpH HyJIeBbIX MOTeHUMasax Hi 
CJBHTOBy!lO — pH 3alaHuv pa3HOCTH NOTeHIMasIOB U HYJIEBOM MeXaHH4eCKOM BO3qecCTBHH. Ucnomb30BaHbl COOTBET- 
CTBYIOIUHe TuMOTesbI OO H3rHOe HM CABHTe, MOCTPOCHEI Be CHCTeMbI OOBIKHOBCHHBIX WHPPepeHUMaIbHbIX ypaBHeHHi 
TpaHW4HBIX YCHOBHH, KOTOpble PeWIalOTCA aHAMIMTH4eCKH 6€3 MCMOb3OBAHHA “TAKEJIBIX) KOHCYHO-3JICEMCHTHBIX MAaKeTOB. 
Ja cpaBHeHHA pe3yJIbTATOB HM MOATBEpYKTCHUA aeKBATHOCTH peIOKeCHHOTO MeTOIa NPOBeeCHO KOHCUHO-3JIEMCHTHOe 
MOJesIMpOBaHHe TakHXx 3a/ja4 B CIelHasIH3HpOBaHHoM Makete ACELAN. 3rTo cpaBHeHue noka3asi0, YTO OMHOKa B OMpesfe- 
JICHHM CMeNIeHHH UM SJIeKTpH4eCKOrO NOTeHHMasa Mp MCHOJ30BaHHU 3TOFO NOAXOA, B Cilyuae 3ajaHuA MeXAaHHYeCKHX 
Harpy30K H pa3HOCTH MOTeHUMaIOB, He lipeBbiiuaeT 6 %. PaspaOoTaHHbiii B CTaTbe METO MOX%KET ObITb IPHMeHEH IIpH IIpo- 
eKTHPOBaHHH Ibe303JIEKTPHYeCKUX TeHepaTOPOB HAaKOIWICHHA IHEPIuH B HU3KOUACTOTHOM OOACTH. 


Konouesbie cJIOBa: YCTpolcTBO cOopa 3Hepruny, TIbe309JIEKTPHAeCKHH reHepaTop, MOpHcTaaA KepaMHKa, w3rH6 
TW1aCTHHBI, CBHI HWIacTHHbI, WpuKilaqHaat TCOpuAa 


Baarojapuocrn. ABTOPpBI BbIpaxKaroT OlaroyjapHocTh pedwakKWHH WypHasia WU peuweH3seHTaM 3a BHUMAaTeJIBHOe 
OTHOIMICHHEe K CTAaTbe. 
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Introduction. Piezoelectric generators (PEG) are used to convert mechanical energy into electrical energy, followed 
by its accumulation. One of the application areas of PEG is the creation of low-power autonomous renewable sources of 
electric energy. The working element of the PEG is a piezoceramic element of a certain shape. The shape and type of 
deformation of this element determine the piezomodule, which characterizes the conversion of mechanical deformation 
energy into electrical energy. Thus, piezomodule d33 is associated with tension-and-compression along the axis of 
polarization, d3; — with the same deformation in the transverse direction to this axis, d:; — with shear. The use of porous 
ceramics makes it possible to create more efficient PEG. This is due to the fact that the elastic modules of porous ceramics 
decrease significantly stronger with increasing porosity than piezomodules. Thus, under the same mechanical load, the 
deformation amplitude of porous ceramics is greater; therefore, the output electric potential is also greater. 

PEG calculation can be performed by the finite element method implemented in ANSYS, ACELAN, COMSOL, and 
others packages. For piezoelectric elements, one or two sizes of which are significantly smaller than others (plates, rods), 
applied calculation theories can be constructed based on hypotheses about the distribution of mechanical and electric 
fields. Without the use of “heavy” finite element packages, applied theories make it possible to model various devices 
based on piezoactive materials. Piezoelectric, piezomagnetic and composite piezomagnetoelectric materials are 
considered as such materials. The construction of these theories is based on the acceptance of hypotheses about the 
distribution of mechanical, electric and magnetic fields. These hypotheses are related to the vibration mode of elastic and 
piezoactive elements of PEG. The most common designs are active and semipassive bimorphs based on multilayer plates, 
polarized in thickness with electrodes on the front faces, performing transverse bending vibrations. A number of papers 
are devoted to the study of devices with shear deformation of piezoelectric elements. An electric model with piezoelectric 
defining equations of mode djs and a single-degree-of-freedom model were combined to describe the energy collection 
characteristics of a piezoelectric cantilever in a shear mode in operation [1]. The proposed model is used to simulate the 
frequency dependence of the output peak voltage and power. The results show a good agreement with the experiment and 
the finite element calculation in ANSYS. In [2], a piezoelectric shear mode energy converter was developed to use the 
energy of a pressurized water flow. It converts the energy of the flow into electrical energy through piezoelectric 
conversion with vibration of the piezoelectric film. A finite element model has been developed to estimate the generated 
voltage of a piezoelectric film, which is in good agreement with the conducted field experiment. A one-dimensional fully 
coupled beam vibration model based on Timoshenko-type hypotheses, which provides a single common basis for energy 
analysis in shear and bending modes, is presented in [3]. In [4], the effect of the inhomogeneity of the properties of the 
plate under shear and torsional vibrations of its central part was studied. In experimental work [5], a multilayer cylindrical 
piezoelectric shear actuator (MCPSA) operating in shear mode ds, was presented for precision actuation under high 
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mechanical load. The actuator was made of piezoelectric ceramic rings Pb(Zr,Ti)O3 (PZT-51), which were concentrically 
assembled together in an electrically parallel connection with alternately positive and negative polarization in the axial 
direction. In [6], metamaterial of identical elementary cells was created, an artificial prototype of a device with 
characteristic patterned electrodes and piezoceramic subunits arranged in a row was designed and manufactured, which, 
as proven, ideally generated the synthetic shear deformation of the face. At the same excitation voltage, there was an 
increase in the displacement of the shear type by more than an order of magnitude, compared to the previous volumetric 
elements in mode djs. In the static formulation in [7], the field of electromechanical coupling in the shear-bending mode 
for an annular piezoelectric plate was theoretically established. In accordance with the classical theory of elastic plates of 
small bending and piezoelectric defining equations, an analytical solution of the bending deformation of the piezoactuator 
under the action of an electric field and a concentrated or evenly distributed mechanical load was achieved. The 
mechanism of generating bending deformation was explained by axisymmetric shear deformation, which additionally 
caused bending deformation of one piezoelectric plate in the form of a ring. This mechanism differs significantly from 
the mechanism of piezoelectric bimorphic or unimorphic drives, which were previously reported. The design of the 
annular piezoactuator has been optimized. In [8], a one-dimensional model was used to construct a sensor response 
function based on shear resonators (quartz cuts) of a bulk acoustic wave, which are promising for in-line measurements 
of fluid viscosity, e.g., in industrial processes. In [9], using the finite element method, a piezoelectric flying height control 
converter was investigated using a shear model deformation. In [10], the theory of a functionally graded plate 
with four-unknown shear deformations was used to express the displacement component. The distribution of the electric 
potential was a linear function of thickness. The plate was under mechanical load and electrical voltage. The basic 
equations and boundary conditions were derived using the virtual work principle. The analysis of stresses and 
deformations from the design parameters was performed. The electromechanical analysis of the stability loss of a 
piezoelectric nanoplate under shear using a modified theory of paired stresses with different boundary conditions was 
studied in [11]. To take into account the electrical effects, an external electrical voltage was applied to the piezoelectric 
nanoplate. A simplified theory of the first-order shear deformation was used. The basic differential equations were 
obtained using the Hamilton principle and nonlinear Von Karman deformations. Finally, the results showed that the effect 
of external electrical voltage on the critical shear load arising on the piezoelectric nanoplate was insignificant. In [12], 
using a combination of two classical approaches to modeling the nonlinear behavior of piezoelectric materials, a 
piezoelectric shear drive for an atomic-force microscope was investigated. Specifically, the novelty of the proposed 
method was in the fact that it combined two sources of nonlinearity of the field-dependent model from Miiller and Zhang 
with the frequency-dependent model from Damjanovic. The numerical results obtained using the finite element method 
(FEM) were compared to the experiment. 

Vibrations in which there is shear in addition to bending are less studied in the scientific literature, 1.e., piezomodule 
d\s, 1s the “working” one, whose value decreases with increasing porosity, but to a lesser extent than elastic modules. The 
latter circumstance makes it possible to build an efficient energy conversion device. Therefore, the development of an 
applied theory of PEG calculation using porous piezoceramics based on simplified models without “heavy” finite element 
packages seems to be a highly relevant task. The objective of this work was to construct an applied calculation method 
for steady-state transverse vibrations in the low-frequency region of a porous piezoceramic plate characterized by both 
shear and bending. 

Materials and Methods. The PEG under study was a piezoceramic plate (length /, thickness /), polarized in thickness, 
cantilevered on the right side, the left side was attached to an inertial mass that performed vertical vibrations and was 
fixed in the horizontal direction. The electrodes were located on the sides of the plate; therefore, with a potential difference 
on them and no mechanical load, the principal deformation in the low-frequency region was shear. Vibrations whose 
frequency was less than the frequency of the first resonance were considered. 

Mathematical formulation of the problem 

The mathematical formulation of the problem is described by a system of differential equations [13] and the 
corresponding boundary conditions. 

Pptitagpju-V-o=f; V-D=0 
o=cj--(€+Byé)—e; -E E=-Vo (1) 
D+caD=e;--(e+Gué)+95-E =(Vu+Vu")/2 


When considering porous ceramics of connectivity 3-0, equations (1) use effective physical constants determined 
through the ACELAN-COMPOS package [14]. These effective properties based on representative volumes (Fig. 1) were 
obtained in [15]. They are presented in Table 1. 
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Fig. 1. Representative volumes in ACELAN-COMPOS package [14] 


Table 1 
Effective properties of porous ceramics 

% of porosity 0 10 20 30 40 50 60 70 80 
P, kg/m? 7,500 6,750 6,000 5,250 4,500 3,750 3,000 2,250 1,500 
ee 10 Nin? 13.90 11.56 9.25 6.85 5.05 3.34 2.07 1.26 0.68 
e, 10”, Nim 7.78 6.15 4.66 3.14 2.10 1.16 0.62 0.28 0.13 
cn , 10%, Nin? 7.43 5.82 4.25 2.82 1.70 1.06 0.52 0.24 0.10 
ee 10" Nia 11.50 9.53 7.23 5.42 3.91 2.72 1.63 0.91 0.47 
chet 10! N/m? 2.56 2.23 1.83 1.44 1.10 0.74 0.44 0.23 0.10 
eo Clin 15.10 13.38 11.37 9.59 7.68 5.93 3.93 2.30 1.25 
es", C/m? 5.20 | -4.23 —3.14 | -2.07 —1.32 —0.75 0.43 —0.21 —0.10 
a Cm? 12.70 10.96 8.96 6.91 5.00 3.30 1.95 1.00 0.44 
Ki, / 86 730 663 582 509 439 349 263 191 122 

ic SEM 635 567 492 413 345 270 199 130 75 


According to Table 1, the dependences of effective elastic modules and piezomodules on the percentage of porosity 
[15] are plotted. They are shown in Figure 2. 
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Fig. 2. Dependences of values on percentage of porosity [15]: a — elastic models; b — piezomodules [15] 
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In accordance with the results shown in Figure 2, elastic modulus £; decreases significantly faster than 


piezomodule ds. 


Building an applied theory 
The developed calculation method consists of solving two problems: in the first, bending under the effect of a 


mechanical load with zero potential difference is considered; in the second, the shear caused by a potential difference in 
the absence of a mechanical load is modeled. In both cases, the absence of charges on the front faces of the plate is taken 
into account. Under the action of a mechanical load and an electrical potential difference, the results of the two tasks add 


up due to the linearity of the task. 
The action of the mechanical load, the potentials on the electrodes are zero. Hypotheses of the Kirchhoff-Love type 


are accepted regarding the equality of normal stresses to zero, and for displacements: 
U3 =UZ> (x), uw --(fuz (*)}2 (2) 
dx 
The distribution of electric potential over the thickness is assumed to be quadratic and symmetrical: 
2 
=D (x)e(22-1]0 +, (x)e(14 22) +, os 4} (3) 


To take into account the boundary conditions at the ends of the plate (x = 0, J), an expression for the transverse force 


is obtained: 


8033233 + 8e33 8 


C13€33 ; re 
1 (cugas +e) 803333 + 8e33 8 
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1 1 8¢33233 +8e33 4 d 
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2 2 
€33 (cisessh — €33€31h ) 


C13 C13 
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+ Cl 
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] eis (cixes3h? ~¢33e31h* )h a 
2 7 UZL2 (x) : 
24 C33 233 + €33 dx 


Taking into account the equality of the normal component of the electric induction vector on the front faces (z = +h) 
to zero, the equations for the unknown deflection UZ2(x) and distribution of the electric potential ®2(x) have the form: 


833233 +8e33 8 


C13€33 7 © 
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1 1 8033233 + 8e33 = 3 ag d : ®, (x) n 
2 2 (casgas +03 )h? h dx 
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The potential difference is set, the mechanical load is zero. Independence of lateral displacement from thickness, 
equality of the longitudinal displacement to zero, and the quadratic distribution of the electric potential over the thickness 


are assumed (3): 
U3 = UZ) (x), m1 =0. (7) 


Expression for the transverse force: 
d d 
Oi =-eu{ Luz, (:) |t-a[ er (3) (8) 


Taking into account the equality of the normal component of the electric induction vector on the front faces (z = +h) 
to zero, the equations for the unknown deflection UZ,(x) and distribution of the electric potential B2(x) have the form: 


2 2 
-ea| Suz o)r-es| So. 9} ~ p(x) -W*phUZ2(x) =0, (9) 
dx dx 


d° d° 
eis a2 (x) el 52 Dr (j=. (10) 


Research Results. The results of calculations based on the proposed applied theories were compared to the 
calculations of piezoelectric element vibrations (/ = 0.1 m, 4 = 0.01 m) at a frequency equal to 100 s" by the finite element 


method in ACELAN [16]. 
In the first problem, defined by equations (5) and (6), when specifying a uniformly distributed load p(x) = 1000 Pa-m 


and with boundary conditions: 
“uz (0) = 0, O1 |r=0= 0, D2 (0) =D, (/) = 0, UZ, (/) = 0, “uz: (/) = 0, (11) 


the following results were obtained, shown in Figures 3, 4. 
The calculations showed that the error in determining the vertical displacement was 5.8%, and for the horizontal 


displacement, it was 1.2%. 
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a) b) 
Fig. 3. Vertical displacement in ACELAN in the first problem: a — distribution; b — graph on the upper boundary 
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Fig. 4. Horizontal displacement in ACELAN in the first problem: a — distribution; b — graph on the upper boundary 


In the second problem, defined by equations (8) and (9), when setting zero load p(x) =0, potential difference 
Vo = 100 W and with boundary conditions: 
O} |x=0= 0, Bz (0) = 0, ®2 (7) =Vo, UZ2 (1) = 0, (12) 


the following results were obtained, shown in Figures 5, 6. 
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Fig. 5. Vertical displacement in ACELAN in the second problem: a — distribution; b — graph on the upper boundary 
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Fig. 6. Electric potential in ACELAN in the second problem: a — distribution; 
b — dependence on the longitudinal coordinate in the middle of the thickness 


The calculations showed that the error in determining the vertical displacement was 0.8%, and for the electric 
potential, it was less than 1%. It should be noted that the values of the horizontal displacement calculated in ACELAN 
turned out to be three orders of magnitude less than the maximum vertical displacement, which indicated the adequacy 
of the hypothesis (7). 

When setting mechanical loads and potential differences, the error of the proposed method turned out to be about 6 % 
for the displacement and electrical potential components. 

Discussion and Conclusion. As noted in the cited literature, the simultaneous use of bending and shear of a piezoelectric 
element can significantly increase its efficiency. In addition, the use of porous ceramics, due to different dependences of elastic 
modules and piezomodules on the percentage of porosity, also improves the output characteristics of PEG. 

In this paper, due to the linear formulation in the theory of electroelasticity, it was possible to develop an applied 
theory for calculating bending-shear vibrations of a piezoelectric element in the low-frequency region, which consisted 
of solving two problems: in the first, mechanical loads operated at zero potentials, and in the second, on the contrary, 
there were zero mechanical loads, and a potential difference was set. Based on various hypotheses about the distribution 
of mechanical and electric fields, two boundary value problems for systems of ordinary differential equations were 
obtained, which were solved analytically. The results of the calculation of displacements and electric potential were 
compared using the proposed method and the FEM implemented in the ACELAN package. These calculations validated 
the applicability of the proposed method, the error for which in calculating the above characteristics was 6%. This 
accuracy is sufficient for engineering calculations; therefore, the proposed method can be applied in the design of 
piezoelectric devices, including the collection and storage of energy. Further development of this applied theory will be 
aimed at covering a wider frequency range, including the first flexural shear resonance. 
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Abstract 

Introduction. Published studies on the rigidity of consoles under load focus on the issues of their deformation and 
destruction. Calculations of the inertia moment, fundamentally important characteristic of the strength of the rod, are 
described. However, the problem of significant time consumption for such calculations has not been solved. The presented 
study meets the lack. The objective of the work is to describe a new rapid method for analytical calculation of the shear 
stress distribution in the section of the console corresponding to the action of an external applied force. For the first time, 
tangential stresses are considered, and examples of calculating the inertia moment for two non-standard sections of the 
console are given in this context. 

Materials and Methods. To develop a new method, the console was presented as a pack of plates oriented parallel to the 
vector of external force. The source calculations were based on the scheme of a console beam with a dedicated plate. The 
deformation of the rod elements was modeled taking into account the effect of a uniform shear stress field in the plate 
section. To validate the simplified calculation of the inertia moment of the sections, schemes of a square, ellipse, triangle, 
hexagon, six-pointed star, and a figured cross were used. Analytical and mathematical research methods were applied, 
specifically, the Huygens—Steiner theorem. 

Results. A rapid valid method for calculating the inertia moment of the cross section of the console under loading has 
been developed. Its difference is the rejection of calculations for each section, taking into account the shape and other 
features. For any shape of the section, the beam is represented as a bundle of infinitely thin plates, their inertia moments 
are integrated, and a well-known solution for deflection of a thin plate is used. The method allows us to unambiguously 
show the distribution of tangential stresses at the end of the console, providing a given deflection, and tangential stresses 
are used for such solutions for the first time. Their profiles are obtained depending on the direction of the external applied 
force. Formulas for the inertia moments of complex sections — a six-pointed star and a figured cross — are derived for 
the first time. Each section is correlated with the stress distribution curve and its maximum value. This data is visualized 
in the form of diagrams. It is found that the inertia moment and the rigidity of the console do not change when the external 
applied force is rotated by 30° for a star-shaped section and by 45° for a square and a figured cross. In general, the tangent 
field depends on the geometry and on the orientation of the section relative to the external applied force. 

Discussion and Conclusion. The proposed simplified approach to calculating the inertia moment of the cross sections of 
the consoles makes it possible to uniquely determine the field of tangential stresses at the end, which provides the 
appropriate value of the external applied force for the given deflection. Engineers and mechanics can use the results of 
the presented work in the calculations and modeling of deformation of rod structural elements. 


Keywords: rod deformation, inertia moment of a flat figure, inertia moment of complex sections, elastic deflection of the 
console, shear stress distribution 
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AHHOTalna 

Beedenue. OrryOmukoBaHHble UCCIeOBaHMA 2KECTKOCTH KOHCONeM MOA Harpy3Koi (OKyCUpylOTCA Ha BOTIpocax Ux Je- 
cbopMallun HM pa3spyuieHua. ONMcaHbl pacueTbI MOMeHTa HHEpuUWH — UpHHUWMMaIbHO BaxKHOM XapakKTepHCTHKH Ipoy- 
HOCTH CTepxKHA. OWHaKO He pellieHa mpoOsIema 3Ha4HTeIbHBIX 3aTpaT BPCMeHH AIA TaKHX BLIEUCIeHHH. IIpeactaBsen- 
Hoe MCCIeAOBaHHe BOCHONHAeT WaHHbIM mpoben. Liem padoTbl — onucaHue HOBOrO ObICTporo MeTOJa aHaIMTHYeCKOTO 
pacueta paciipeyesIeHua HallpmKeHHA CABHTa B CCHCHHM KOHCOJIM, COOTBETCTBYIOM[ero AeHCTBHIO BHeIUHeH MpHy10x%KeH- 
HOM cnbt. BrepBbie B TAKOM KOHTeKCTe pacCMaTpHBaIOTCA KaCaTeJIbHble Hallps.KCHHA H IPHBOAATCA IIpHMepbl pacueta 
MOMEHTa HHepwHH Jit TByX HeCTaHJapTHbIx CeyveHuii KOHCOJIN. 

Mamepuaaoi u memoodot. na co3qaHud HOBOTO MeTOa KOHCOJIb MIpeACTaBMJIM Kak Madky WWlaCTHHOK, OPHeHTHpOBaH- 
HBIX TlapaJIIeIbHO BeKTOpy BHeWHel cub. UcxoWHble pacueTbI CTPOWJIM MO CX€Me KOHCOJIBHOM OasKH C BbIAeeHHOM 
TIacTHHKON. JecbopMativ1o CTep2KHEBBIX IJIEMCHTOB MOJeIHpOBasIH C YAeTOM JelCTBHA OAHOPOAHOTO MOA HalIpsKe- 
HHA CBHTa B CCUeHHH TWIacTHHKH. Jd OOOCHOBaHHA YIPOWIeHHOTO pacueTa MOMEHTa HHeEpUHH CeueHHi 3a ,eHCTBOBaIN 
CX€MBI KBaj[pata, 9JWIMTICa, TpeyroJIbHUKa, WeCTHYTOJIbBHUKa, WICCTHKOHCYHON 3Be31bI H (urypHoro Kpecta. Mcnomb30- 
BaJIM AHaJIMTH4eCKHe HW MATeEMATHYeCKHEe METOAbI HCCIeOBaHHA, B YaCTHOCTH Teopemy I toMrenca—lLTetinepa. 
Pe3yivmamot ucciedoeanua. Co3yaH ObIcTpbIii YHHBepcasIbHBIM MeTO], BBIYHCICHHH MOMeHTa HHeEpuHH MonepedHoro 
Ce4YeHHA KOHCONM HOA Harpy3kou. Ero ormmune — OTKa3 OT pacueTOB JIA KaxKOTO CeYeHHA C y4eTOM (POpMBI H Apyrux 
ocoOeHHoctett. [pu mobo dopme ceyeHua Oalka WpeycTaBsaeTcA Kak Nauka OeCKOHeYHO TOHKHX IWIaCTHHOK, MO- 
MCHTbI HX HHEPLMH MHTerpupylorca, M HCMONb3yeTCA U3BECTHOe pellieHHe AJId WporuOa TOHKOH MWacTHHKU. Metoy 1103- 
BOJIACT OHO3HAYHO MOKa3aTb paciipeseseHve KacaTesIbHbIX HallpsKeHHH Ha TOPLe KOHCOMH, OOeCHeYMBaIOLNHX 3aaH- 
HbIi MIporHo, MmpH4em BIepBBIe JIA TAKHX PeLWICHHM HCIONb3Y!IOTCA KacaTeJIbHble HalilpsKennaA. Ilomyyenbl ux npopuru 
B 3ABHCHMOCTH OT HallpaBsIeHHA BHELIHeM IPWIOKCHHON cHIbI. BriepBble BEIBECHbI POPMYVIbI JIT MOMCHTOB HHepunH 
CJIOXKHBIX CeUeHHH — WICCTHKOHCYHON 3Be3AbI U (urypHoro Kpecta. Kaxxyoe ceyeHue COOTHECeHO C KpHBOM pacripeye- 
JICHHA HalipsxKeHHA HM ETO MAKCHMAJIBHBIM 3HadeHHeM. ITH JaHHble BU3YaJIM3HpPOBaHbl B BUe WMarpaMM. YcTaHoBJIeHo, 
4TO MOMCHT HHEPIHH WM 2KCCTKOCTL KOHCOJIM HE MeCHAIOTCA IPH MOBOpOTe BHEMIHel IpHIORKeHHOM cnt Ha 30° pA 
CeYeHHA B Be 3Be3bI H Ha 45° — yA KBay{pata u (purypHoro Kpecta. B oOujem cyryyae Moyle KacaTeJIbHbIX 3aBHCHT 
OT rCOMeTPHYeCKO (POPMBI HU OT OPHEHTAL[HH CeYCHHA OTHOCHTEJIBHO BHELIHEM MpHI02%KCHHOM CHIBI. 

O6cyotcoenue u 3akmo4uenue. IpeqioxKeHHbI yIpOWeHHbIi MOAXO K pacueTy MOMeHTa HHEpLM NOMepedHBIx ceue- 
HHH KOHCOJIeH WaeT BO3MOXKHOCTS OJHO3HAYHO OMIPeeIMTb NWO KaCaTCIbHBIX HallpsKeHUM Ha TOple, OOecreuHBalO- 
wee TIpH 3aaHHOM Tpornbe cooTBeTcTByIOlce 3HadeHHe BHeELMIHeH MIpH1O%*KeHHON cub. WHoKeHepbI MH MCXAHMKH MOTyT 
MCIOJIb30BaTb Pe3yIbTAThI MpeACTAaBIeHHOM paOoTb! Up pacueTax H MOeIMpoBaHHH JePopMalHM CTep2KHEBBIX 9JIe- 
MCHTOB KOHCTpyKUHit. 


Karouesble C10Ba: JepopMalluaA CTepxKHA, MOMCHT HHeEPIHH TWIOCKOM PUrypbl, MOMeHT HHEPLUMH CiIOXKHBIX CeYeHHH, 
yiipyrHii mporm6 KOHCOJIM, paciipeyemeHue KacaTeJIbHbIX HallpskKeHHii 
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Introduction. Numerous building structures contain elements in the form of rods, which undergo elastic deformations 
during manufacture or operation [1]. The bending strength of a rod or beam determines the bearing capacity of the 
structure [2]. The ability of a beam to elastic deformation is characterized by stiffness, defined as the ratio of load P to 
the elastic deflection of the beam 4. : mx = P/A- [3]. As a rule, under laboratory conditions, stiffness is checked on a 
console beam. One end of it is embedded in a rigid base, and an external force, directed perpendicular to the axis of the 
beam, acts on the other one [4]. Modeling and calculations of deformation and fracture characteristics of rods, as a rule, 
are associated with solving the problem of deflection of a console beam, or console, under the action of an external applied 
force [5]. At the same time, there are no publications on simple and valid methods for determining the inertia moments 
of complex sections relative to the action of an external applied force. The solution to this problem is described in the 
presented article. 

This work is aimed at the creation of a valid, rapid method for calculating the inertia moments of complex console 
sections under the action of an external applied force. The new approach provides analytical determination of the shear 
stress distribution in the section corresponding to the action of an external applied force. It should be noted that earlier, 
tangential stresses were not taken into account in such calculations. In addition, for the first time, examples of calculating 
the moment of inertia for complex figured sections are given. 

Materials and Methods. In a number of works on the resistance of materials, e.g., in [6], a universal formula is given 
for calculating the elastic deflection of console Ac. According to this formula, the rigidity of the console is: 

P/.,, =3EI,/D, (1) 
where E — Young's modulus; Z — length of the console; /, — inertia moment of the cross-section of the beam relative 
to x-axis, passing through the center of gravity of the section perpendicular to the applied force P. 

It follows from equation (1) that a fundamentally important characteristic of the console is the inertia moment of 
section J,, whose value depends on the geometry of the cross-section of the beam and the direction of x-axis [7]. It 
should be emphasized that in equation (1), inertia moment J; refers to x-axis, which is perpendicular to the direction of 
the external applied force P. Specifically, the inertia moment of rectangular section axb relative to the axis of symmetry 
x is equal to [8]: 

I, = ab? /12. (2) 

Here, a — thickness of the console, b — its width. Force P is directed parallel to side b of the rectangle. 

Substituting (2) into (1), we obtain the well-known equation for the deflection of a rectangular console [9]: 


3 
oO 
Ea\ b 


The cross sections of console beams, or consoles, are different. Figure | shows a simple example of a square section 
console under the action of external force P, directed along the diagonal of a square. 


ve 


Fig. 1. Diagram of a console beam with a dedicated plate 


The inertia moment of the section relative to x-axis is called the sum, or integral, of the products of elemental areas 
ds = dxdy by the squares of distances y of the areas to x-axis: J, = {J y*dxdy. [10]. The integrand function is virtually the 


inertia moment of the elemental area dxdy relative to x-axis. 

The console can be represented as a pack of extremely thin plates with thickness dx and length LZ, oriented parallel 
to the force vector P. Under the action of P, all plates bend by the same amount A,. With a given orientation of the 
console section, a separate plate is not affected by deformation of the rest of the volume. Then, the inertia moment of 
the section of the console as a whole will be determined by the integral sum of the inertia moments of the sections of 
all the plates in the pack. 

The projection of the plate onto the plane of the console cross section is a rectangular strip with thickness dx and half- 
length h (Fig. 1). The inertia moment of the section of a separate plate can be considered as the inertia moment of the 
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console of rectangular section 2hxdx. By definition, the equation of type (2) is applicable to each plate in a pack, where 
a= dx and b= 2h. According to this expression, the inertia moment of the strip is dl. = 2/°(x)dx/3. Thus, the inertia 
moment of the console section can be determined through integrating the inertia moments of elementary strips rather than 
elemental areas: 


i= = [7 AGo at, (3) 


where x varies from A to B. 

The condition for the orientation of the plane of the plates parallel to the vector of the external applied force is 
important, since it provides a unique association of the elastic deflection of console 4. with the distribution of tangential 
stresses in the cross section of the console. All the plates in the pack bend by the same value ).. according to (1), 
dl, = 2h?(x)dx/3. This means that elementary force dP = 3/.EdI,/ L? = 21-Eh*(x)dx / L? is required for a plate with thickness 
dx. This force corresponds to the action of a uniform shear stress field in the plate section ds(x) = 2h(x)dx: 

t=GP/dss1,Eh (iL. (4) 


It can be seen from (4) that the value of voltage A, in the coordinate system xy does not depend on coordinate y. 

Equation (4) is convenient to be used when modeling the deformation of rod structural elements. 

Integral (3) determines the inertia moment of section /, relative to x-axis, passing through the center of gravity of the 
section. In the case of asymmetric and complex sections, it is convenient to first find the inertia moment of the section or 
part of the section relative to the axis that does not pass through the center of gravity of the section. Then, you need to 
move on to the inertia moment of section /, relative to the axis that passes through the center of gravity of the section. It 
is known that the inertia moment of the section repeats the properties of the inertia moment of a solid and obeys the 
Huygens-Steiner theorem [11]. The inertia moment of section J, relative to arbitrary x-axis is equal to the sum of the 
inertia moment of this section /, relative to the axis passing through the center of gravity of the section parallel to x-axis, 
and the product of the cross-sectional area S by the square of distance a between the axes: J, = J, + a*S. Therefore, in the 
general case, we can write: 


_ 258 3 De. 2 
I= 3], h(xdx+a’S =I, +a’S. (5) 


If x-axis passes through the center of gravity of the section, then distance a = 0 and equation (5) turns into (3). 

Research Results. To create a simple, fast, valid calculation method, we abandon calculations for each section, taking 
into account its shape and other features. This approach was implemented for the first time in the framework of this 
research. No matter how complex the cross-section is, it is quite sufficient to use a well-known solution for deflecting a 
thin plate, present the beam as a pack of infinitely thin plates, and integrate their inertia moments. In addition, the method 
allows us to unambiguously show the distribution of tangential stresses at the end of the console, providing a given 
deflection. It should be emphasized that tangential stresses are considered in this context for the first time. 

Short-Cut Calculation of the Inertia Moment of Simple Sections. To validate the proposed method, we consider 
the known sections of simple geometry. Next, instead of the expression “inertia moment of the cross section of the 
console”, we use the term “inertia moment”. We assume that the external applied force is always directed perpendicular 
to x-axis, relative to which the inertia moment of the section is determined. 

Square. Equation (2) is obtained under the condition that the vector of force applied to the end of the console is 
perpendicular to side a of the rectangle. For b = a, we obtain the inertia moment of the square section J, = a4/12. 

Using the proposed method, we find the inertia moment of the square relative to x-axis, which is parallel not to the 
side, but to the diagonal of the square (Fig. 2 a). 


Bd t, MPa 
30 
20 
% 

a 10 

Ob = 

=3 —2 =1 0 1 2 x, mm 

a) b) 


Fig. 2. Cross-section diagram in the form of a square: 
a — cross-section; b — voltage distribution t along x-axis 


Advanced Engineering Research (Rostov-on-Don). 2024;24(2):159-169. eISSN 2687-1653 


The half-length of the strip in Figure | is equal to A = x. From equation (3), we find: 
2 a/3 


c== aes 

3 J-a/i3 12 

It can be seen that turning x-axis by 45° does not change the inertia moment of the square section. Consequently, when 
the external applied force is rotated by 45°, the rigidity of the console does not change either (1). 

Substituting x in place of h (4), we obtain a shear stress distribution in the section of the square console corresponding 
to deflection A: 


I 


Ux) =A,°E/ DP. 

Figure 2 b shows distribution t (x) along the diagonal of the square at A. = 2 mm and a = 5 mm. According to (1), the 
action of the shear stress t (x) on the end of the steel console (FE = 200 GPa) with length Z = 50 mm corresponds to the 
action of an external applied force P = 0.25EA.a‘/L? = 500 N. Maximum voltage tmax = 40 MPa is observed on the vertical 
diagonal of the square. Under the deviating from the diagonal, voltage t decreases sharply to zero. As a result, a sharp 
peak is formed in the system T (x). 

Voltage t depends only on variable x; therefore, according to the graph in Figure 2 5, it is possible to determine value 
t at any point of the square. 

Obviously, when the force is oriented perpendicular to the sides of the square, t = P/a” = 20 MPa (Fig. 2 5, dotted 
line). As can be seen, stress distribution t(x) in the cross section corresponding to external force P significantly depends 
on its direction. 

Ellipse. Figure 3 a shows a diagram of an ellipse with semi-axes a and b. The origin of coordinates is in the center 


of the ellipse. 


t, MPa 


60 


40 


qd) b) 


Fig. 3. Section diagram in the form of ellipse with semi-axes a and b: 
a — ellipse; b — distribution t along semi-axis a. 1 —a=b, 2— b =2a,a=2.5 mm 


It follows from the canonical equation of ellipse [12] that the strip highlighted in Figure 3 has half-length 
h=b(a@ —x’)'/a. Substituting this value in (3) and integrating from —a to +a, we obtain the inertia moment of the 


elliptical section: 


é 2b° 7 3/2 b 

I,=-| Wdx= = [ 2? -x? | ice. (6) 
3 J-a 3a” Ja 4 

At a=b =r, we derive the inertia moment of the circular section: 27“/4, where r — radius of the circle. 


Substituting value / corresponding to the ellipse in (4), we obtain the following distribution of shear stress in the 


section of the elliptical console: 
U(x) =A,b°(1 —x? /a*)E/L. 

Figure 3 b shows distribution t(x) along semi-axis a = 2.5 mm at A, =2 mm, F = 200 GPa and L=50 mm. When 
comparing it to distribution t(x)in Figure 2 b, for a square, a significant influence of the geometry of the console section 
on the stress distribution is obvious. In the case of the elliptical section, there is no sharp peak on curve t(x). Maximum 
voltage Tmax 1s observed along semi-axis b. Taking into account the requirement i. = const, it can be said about a rapid 
growth of voltage Tmax with an increase in ratio b/a. When half-axis b is doubled, tnax increases by 4 times (Fig. 3 5). 
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Triangle. Consider the section in the form of an isosceles triangle (Fig. 4). x-axis is directed along the altitude of the 
triangle (Fig. 4 a). At distance x from the base of the triangle, the half-length of the strip is equal to = b (a —x)/2a. 


ah 
a 
x 
a > b , 
) b) 


~< 
a 


Fig. 4. Section diagram in the form of isosceles triangle: 
a — x-axis, force is directed parallel to the base; 
b — force is directed perpendicular to the base 


Integrating expression (4) over x from 0 to a determines the following value of the inertia moment: 
an bs a, a 
oa [ (a-x}ide= (7) 
Consider the case when x-axis passes through the center of gravity and is parallel to the base of the triangle (Fig. 4 5). 
The strip at distance x from the left corner of the triangle has half-length = xV4a? /b? -1 /2 (Fig. 4b). The center 
of gravity of the strip is located at distance h from the base of the triangle (Fig. 4 b). Therefore, according to the Huygens— 
Steiner theorem for the cross section, the inertia moment of the strip relative to the base of the triangle is: 


3 
di, = 7 Wdx+ 2h de = 2% wa. 
3 3b 


I,= 


Integrating the resulting expression over variable x from —b/2 to b/2 determines the inertia moment of the triangle 
relative to the base: 
1, =a bs 12, (8) 
The center of gravity of the triangle is located at distance a/3 from the base. According to the Huygens—Steiner 
theorem [13], the inertia moment of a triangle relative to its own center of gravity is less than (8) by: 


where S = ab/2 — area of triangle. 
In this manner, the inertia moment of the triangle relative to its own center of gravity is equal to: 

_ab ab _a’b (0) 
12 18 36 

The result obtained corresponds exactly to the tabular value of the inertia moment of the section relative to the axis 


through the center of gravity parallel to the base of the triangle. 
Figure 5 shows the shear stress distributions in the triangular section at A. = 2 mm for cases where the deflection force 


is directed along the base (a) and along the altitude of the triangle (5). 


I, 


t, MPa 


t, MPa 
15 30 


10 20 


10 


Fig. 5. Distributions t in the triangular section: 
a — force is directed along the base; b — force is directed along the altitude of the triangle 
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When the triangle is oriented as in Figure 4 a, voltage t(x) decreases gradually from tmax at the base to 0 at the 
apex (Fig. 5 a). With the orientation of the triangle as in Figure 4 5, distribution t(x) (Fig. 5 5) is similar to the distribution 
for a square section when the force is oriented along the diagonal (Fig. 2 b). 

Regular hexagon. Figure 6 shows two orientation of the hexagon relative to x-axis: parallel (a) and perpendicular (6) 


to its diagonal. 
dx 


b) 


Fig. 6. Regular hexagon: 
a — diagonal is perpendicular to the force vector; 
b — diagonal is parallel to the force vector 


The section of the hexagon in Figure 6 a consists of the following: 

— rectangle with width a and height a3, 

— two isosceles triangles with apex a/2 and base avV3. 

Let us find the inertia moments of the specified parts of the section by (2) for the rectangle and by (7) for triangles. 
For the rectangle in equation (2), b = a3; therefore, its inertia moment is J: = V3a4/4. For triangular parts in equation (7), 
the base is b = aV3, and the altitude is a/2. Consequently, the inertia moment of the hexagon is: 

I, = Ty +19 = V3a* (1/14 41/16) =5V3a4/16. (10) 


If the applied force is directed along the diagonal of the hexagon (Fig. 6 5), then the half-length of the strip is equal 
to a/2 — x/V3. Having made the appropriate substitutions in (4), we obtain: 


2 pav3/2 ; 5 a) as 
I,== an x/ 3 Bdx = re (11) 
Comparing (10) and (11), we make sure that the rotation of the hexagon by 30° does not affect its inertia moment 
relative to x-axis. 
Consider the case of the hexagon orientation as in Figure 6a. At X}=2 mm and £ = 200 GPa in the range from 
—a/2 to +a/2, voltage T(x) is constant and equal to Tmax = 15 MPa (Fig. 7 a). Under the same conditions and orientation of 
the hexagon as in Figure 6 5, distribution t (x) is similar to the distribution for a square section (Fig. 2 b). However, here, 


Tmax = 20 MPa and Tmin = 5 MPa (Fig. 7 b). 


t, MPa t, MPa 
127 16 
oF 12 
67 8 
3} 4 
0.0 1.0 2.0 3.0 x,mm 0.0 1.0 2.0 x, mm 
a) b) 


Fig. 7. Distribution (x) in the section of the regular hexagon: 
a — diagonal is perpendicular to the force vector; 
b — diagonal is parallel to the force vector 


Examples of Simplified Calculation of the Inertia Moment of Complex Sections 

Six-pointed star. Using the proposed method, we calculate the inertia moment of a non-standard section in the form 
of a regular 6-pointed star with side a. To determine the inertia moment relative to the minor diagonal of the star, we use 
the diagram in Figure 8 a. Let us highlight three zones in the diagram: zone I with width a/2, zone II for the rest of the 
half-figure, and the adjacent auxiliary zone III in the form of a triangle. 
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a) b) 


Fig. 8. On calculating the inertia moment of the six-pointed star relative to a — minor diagonal; 
b — major diagonal 


Half-length of the strip in zone I at distance x from the diagonal is 4; = V3(a —x). Using expression (3), for two zones I 
in Figure 6 a, we obtain the inertia moment: 
— wf (a —x)*dx = 15-3 oa 
a/2 16 
The inertia moment of zone II is equal to the inertia moment of a rectangle with height aV3 and width a without the 
inertia moment of triangle III. The half-length of the strip in the rectangle is equal to Ay = a V3/2. From equation (3), we 
find the inertia moment of the rectangle: 


3 
Tian = a. 


4 
The half-length of the strip in the triangle is equal to 4y = xV3. According to (4), the inertia moment of the triangle is: 


v3 
lin = a". 


32 
Hence, the inertia moment of zone II is equal to: 


3B 4 
Tq = Lean — Lin ~ 35 4 : 
The doubled sum of the inertia moments of zones I and II determine the inertia moment of the 6-pointed star: 
11-3, 
ar a a’. 

To determine the inertia moment of the star relative to the major diagonal, we use the diagram in Figure 8 b, where 
two zones are highlighted. The half-length of the strip in zone I is hy = a + x/N3, in zone II — hy = a/2 — x3. From (3), 
the inertia moments of zones I (J = a465V3/96) and II (Jur = a*V3/96) are calculated. It is seen that the doubled sum of 
the inertia moments of zones I and II corresponds exactly to equation (12). Therefore, the inertia moment of the 6-pointed 
star relative to the major and minor diagonals is the same. 

Figure 9 a shows the shear stress distribution in the section of the 6-pointed star corresponding to 1 = 2 mm and the orientation 
of the external applied force according to Figure 8 a. Along the vertical axis of the star, the voltage takes on maximum value 
Tmax = 37.5 MPa. With distance from the axis to the right, t drops quickly to the level t= 9.375 MPa and remains constant in the 
range 2 <x<4mm. Then, t drops to zero. With distance away from the axis to the left, t changes similarly. 


Iq (12) 


t, MPa Tt, MPa 


Fig. 9. Distributions t in the hex-shaped section: a — force is directed along the major diagonal; 
b — force is directed along the minor diagonal 
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When the star rotates by 30°, voltage distribution t changes significantly, acquiring the outlines of a dovetail (Fig. 9 5). 
Two peaks are observed with value Tmax = 28.12 MPa. The distribution pattern is symmetrical. However, with distance from 
the axis of symmetry, the voltage first increases from 12.5 to 28.12 MPa, and then drops to 3.125 MPa. In this regard, two 
peaks are observed on dependence t(x). Further on, the voltage decreases rapidly to zero. 

Figured cross. Consider a solution for the inertia moment of a non-standard section in the form of a figured 
cross (Fig. 10), each side of which is the fourth part of the circle of radius R. The distance between the apices of the figure 
is RV2. We first find the inertia moment of the cross relative to the major diagonal (Fig. 10 a). 


a) 


Fig. 10. Figured cross: a — force is directed along the diagonal of the cross; 
b — force is directed at an angle of 45° to the diagonal of the cross 


Half-length of strip / in Figure 10 a is equal to: 


h=R-.JR?-—(R-x)?. 


Taking into account the symmetry and using equation (3), we obtain the following value of the inertia moment of the 
section of the figured cross: 


1.=4["[R Jx(2R D)] d= [4-5/4] 4 (13) 


3 

We determine the inertia moment of the cross when it is rotated by 45° relative to x-axis (Fig. 10 5). 

Figured cross fits into the circle of radius R. Figure 10 5 shows that there are four figures in the form of an oval with 
sharp corners around the cross. To determine the inertia moment of the cross, it is sufficient to subtract the inertia moments 
of the four ovals from the inertia moment of the circle. 

It follows from (6) that the inertia moment of the circle is equal to: 

To = 1R*/4. 
We determine the inertia moments of the ovals relative to the center of gravity of the cross. 
The inertia moment of the oval on the right is equal to the inertia moment of the oval on the left. Half-height of the 


h, = JR? -(R-x)’. 


We take into account equation (3), as well as the symmetry of these parts and their location. By integrating from 0 to 
R—R/\2, the following value of the inertia moment of these two parts is obtained: 


3 
ie Vi fax RF ] v= 2 2-2 (14) 
3 Jo 3,4 

The centers of gravity of the ovals above and below the circle are located at distance a = R/V2 from the center of 
gravity of the whole figure. The area of these two parts is equal to S = R?(x — 2). 
By (4), we calculate the moment of inertia of this pair: 
Ixy =a?S = R*[n/2-]]. 


strip at the oval on the left is: 


Value h in (3), according to Figure 10 5, (oval at the bottom of the circle) is equal to h = 2a? —x” —a. Therefore, 
for a pair of ovals at the top and bottom of the circle, the inertia moment relative to their own centers of gravity can be 


R/Y2 . 
i= =5Rf [i-(e7Ry 1/3 | dx = R*|[3n/4-7/3]. (15) 


0 


recorded as: 
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We determine the resulting inertia moment of the figured cross. To do this, we subtract the inertia moments of the 

four ovals from inertia moment (13) of the circle: 
Tq Ly, —Le2 — 2153 =[4—500/4] R*. 

Comparing the result obtained to result (13), we can see that the inertia moment of this cruciform section of the console 
does not change when x-axis is rotated by 45°. 

Figure 11 a shows the distribution of shear stress in the cross section of the figured cross, corresponding to A = 2 mm, 
and the action of an external applied force along the axis of the cross (Fig. 10 a). A sharp voltage peak is observed on the 
vertical axis of the cross, where Tmax = 80.0 MPa. 


t, MPa t, MPa 
R=2mm 

60 16 

12 
40 

8 

20 4 

0 0 

= 0 2 x,mm 3 2 -1 0 1 2x,mm 
a) b) 


Fig. 11. Distributions t in the section of the figured cross: 
a — force is directed along the diagonal of the cross, 
b — force is directed at an angle of 45° to the diagonal of the cross 


When the cross is rotated relative to the applied force by 45°, voltage distribution t changes significantly. As the axis 
of symmetry deviates, voltage first gradually increases from t= 13.73 MPa to t max = 20.32 MPa. Then, voltage t drops 
sharply to zero. Therefore, two symmetrical peaks (Fig. 11 5) are observed in distribution t at a distance of 2 mm from 
the axis of symmetry. 

Discussion and Conclusion. The proposed simplified calculation method provides for the rapid determination of the 
inertia moments of complex cross sections of the console. In this case, the shear stress field in the sample section 
corresponding to the action of an external applied force is uniquely determined. In addition, it is shown that the stress 
distribution in the section qualitatively and quantitatively depends on the orientation of the section relative to the direction 
of the external applied force. 

To validate the method, the inertia moments were calculated not only for known sections of simple geometry (which 
showed the absolute identity of the calculated and published results in the literature), but also for two new complex 
sections in the form of a regular six-pointed star and a figured cross. It is shown that the rigidity of the console does not 
change when an external force applied perpendicular to the axis of symmetry is rotated by 30° for a section in the form 
of a 6-pointed star, and by 45° — for a square and a figured cross. The method and the solutions obtained can be used by 
engineers and mechanics in modeling and calculating strength and stiffness of rod structural elements. 
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Abstract 


Introduction. The development of the polar areas of the World Ocean and the need to solve various problems 
associated with a large number of freezing inland water bodies issue new challenges for science. These challenges 
include the problem of studying the behavior of ice cover when exposed to various types of loads. Of great interest is 
the consideration of problems about the action of a moving load on the ice cover. A moving load simulates the effect 
of moving vehicles on ice. However, in papers devoted to the above problems, cases of load movement along a straight- 
line trajectory are considered. The objective of this research is to develop a method for studying the behavior of ice 
cover under the action of a load moving arbitrarily. 

Materials and Methods. The article proposes a method for solving the problem of the action of a force moving along 
an arbitrary trajectory on the ice cover of a reservoir of finite depth. The problem amounts to solving a system of two 
differential equations. The first of them models the behavior of the ice cover, and it is the equation of vibrations of a 
viscoelastic plate. The second equation simulates the behavior of fluid in a state of potential flow, and it is Laplace's 
equation. To solve the system of differential equations, integral transformations in time, space and variables were used. 
The resulting solution was expressed through an iterated integral, which was calculated using numerical methods. 
Results. The development and implementation of the method resulted in solving the problem of the movement of a 
concentrated force along an ice cover according to an arbitrary law. At the same time, studies were carried out on the 
behavior of displacements and stresses in the ice cover depending on the speed and acceleration of the movement of 
the vertical load, on the depth of the reservoir, and on the viscoelastic properties of ice. In addition, the distribution of 
the velocity vector of fluid particles along the depth of the reservoir was calculated. 

Discussion and Conclusion. The proposed method is very effective for solving problems of moving loads acting on 
the ice cover of a reservoir of finite depth. It provides solving problems about the action of a load moving along an ice 
cover along a complex trajectory. The results obtained can be used to calculate the stress and displacement of the ice 
cover during the laying of ice roads or the construction of airfields on the ice. 


Keywords: infinite ice cover, moving load, arbitrary trajectory, variable speed 
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Pa3pa6oTKa MeTOJa pellleHHA 3afa4n AepopMauHn JeqAHOrO NOKpPOBaA TOA eucTBHeM 
IIPpOW3BOIbHO JBIOKyWeHca Harpy3KH 
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AHHOTauna 

Beedenue. Ocsoenue nouApHbIx pationos Muposoro okeaHa, HEOOXOAMMOCTb pellieHHA pa3ssIM4HBIX 3aja4, CBA3AHHBIX 
C HaMyveM OONbWIOTO YMCA 3aMep3aIOWHX BHYTPCHHHX BOJOCMOB, CTAaBAT Tepe HayKOl HOBbIe mpoOsempl. K ux 
YHCIIy OTHOCHTcA MIpoOIeMa H3YHCHHA NOBeACHHA JeAAHOTO MOKpoBa MOA BO3AeHCTBHeM Ha Hero pa3IM4HOrO Bua 
Harpy30K. bombo HATepec NpeAcTaBsIAeT paccMOTpeHHe 3afay O elCTBHM Ha JIeqAHOM MOKPOB NOABWKHOHM Harpy3Ku. 
TloqBpywxKHad Harpy3ka MOjjeMpyeT AelicTBue Ha JIeq ABIMDKYIMXCA TpaHCOpTHbIx cpeycTB. OyHako B paodoTax, 
TIOCBAL[CHHBIX BbIMIeyKa3aHHbIM 3alayaM, paCCMaTpHBaIOTCA CilyyaH BMWKeHHA Harpy3KH MO UpaAMOMMHeHOH 
Tpaextopuu. Llenbro anno paooTsl ABIAeTCA pa3spaOoTKa MeTOJa HCCIeAOBAaHHA MOBeACHHA JIeAAHOTO MOKpoBa Toy 
TeHCTBHeM Harpy3KH, MepeMellaroulelica IPpOM3BOJIbHEIM OOpa30M. 

Mamepuaavi u memoool. B cratbe ipeqoxKeH MeTOX pelleHHA 3aa4uH O AevCTBHM Ha JIeqAHOM MOKpOB BOoeMa 
KOHeYHON PIYOHHEI ABMOKYWeHCA 10 IPOM3BOJIbHOM TpaeKTOpHH CuIbI. 3ayaya CBOAMTCA K PeIeHHIO CHCTeMBI JByX 
TuddepenunanbHEix ypaBHennn. Ileppoe 13 HUX MOJeMpyeT MOBeeHHe eqAHOTO TOKPOBa H ABJIACTCA ypaBHeHHeM 
KoseOaHHi BA3KOYUpyrou WiacTHHbI. Bropoe — MoyesMpyeT NOBeeHHe 2%KHAKOCTH, HaxOAulelica B COCTOAHHU 
MOTCHUMaIbHOrO TeYeHHA, MH ABIAeTCA ypaBHeHuem Jlanmaca. JIna pewleHua cHcTeMbI AUdpdepenuMabHBIx 
ypaBHeHHM TIPHMeHAIMCb WMHTerpasibHble MpeoOpa30BaHHA MO BpeMeHHOM HW MpOCTpaHCTBeHHbIM TIepeCMeCHHBIM. 
IlomyueHHoe B pe3yIbTaTe pelieHHe BbIpaxkaocb Yepe3 MOBTOPHBIM UHTerpall, AIA BbIYMCIeHHA KOTOPOTO 
IIPHMCHAJIMCb YACICHHBIC METOBI. 

Pe3yismamoet ucciedoeanua. B pe3yibTate peaM3alMn WpesoKeHHOrO MeTOJa HoMyYeHO pelleHve 3aqaun oO 
TBWKCHHH COCpexOTOYCHHOM CHJIbI 10 JIeqAHOMy MHOKpoBy MO MpOu3sBOIbHOMy 3aKoHy. IIpu 3TOM Mpou3BeeHbI 
McCCIeqOBaHHA XapakTepa MOBeCHHA MepeMeICHHH M HallpsKeHHH B JIe{HOM MOKPOBE B 3ABHCHMOCTH OT CKOpOCTH H 
YCKOpeHHA JBMOKCHHA BeEPTHKAIbHOM Harpy3KH, TyYOMHbI BOOeCMa UM BASKOYUPyrux CBOMCTB JIbya. Kpome toro, 
paccunTaHo paciipeyjeeHve BeKTOpa CKOPOCTH YacTHL *XUAKOCTH M0 rryOuHe BoOeMma. 

O6cystcdenue u 3akiio4enue. [penioKeHHbii MeTO ABIAeCTCA BeCbMa 93*PPeKTHBHBIM JIA PeWIeHHA 3aa4 oO 
MOABWOKHBIX Harpy3Kax, J[eHCTBYIOWIHX Ha Je AHO MOKPOB BOJOeMa KOHeYHO! TryOuHbI. OH MO3BOIAeT pellaTb 3aauH 
0 WeHCTBUM Harpy3KH, JBMKyWelica 10 eqAHOMY MOKpoBy M0 cIOX%KHOM TpaexTopun. IosyuenHble pe3syIbTaTbl MOryT 
ObITb HCHOJb3OBaHbI JIA pacueTa HallpMKeHHA HM TepeMeLeHHH JIEMOBOTO MOKpoOBa Mp MpOKaAKe JIEOBbIX JOpor WH 
CTPOHTEJIBCTBE aIPOAPOMOB Ha JIBALY. 


K.ouespie cs10Ba: OecKOHeUHBII Te qHOU TIOKPOB, ABMWKylWaAcd Harpy3Ka, MIpOW3BOJIbHad TpaecKTOpHA, MepeMeHHadA CKOPOCTb 


BuaarogapHocru. ABTop BbIpaxKaeT OaroapHOCTB peueH3eHTaM 3a YKa3aHHble 3aMe4aHHA, KOTOPbIe NOSBOJIMJIN 
TIOBbICHTb KAYCCTBO CTATbH. 


Aisa wutTupospannsa. Tanadypqun A.B. PaspadoTka MeToga pellieHua 3ayqa4uu ePpopMaluu seqaHoro NOKpoBa Toy 
TeVicTBHeM MpOv3BONbHO ABMKylWelica Harpy3ku. Advanced Engineering Research (Rostov-on-Don). 2024;24(2):170-177. 
https://doi.org/10.23947/2687-1653-2024-24-2-170-177 


Introduction. The development of the polar areas of the World Ocean and a large number of freezing inland 
reservoirs drive the need to study the fields of displacement and stresses of the ice cover caused by the action of various 
types of loads. Numerous papers by domestic and foreign scientists are devoted to solving these problems. Previously, 
it has been found that the mechanical properties of ice depend on its temperature and water salinity. Much attention 
was paid to the development of numerical ice models that accurately reflected the interaction of ice and ideal 
incompressible fluid. In [1, 2], the smoothed particle hydrodynamics was used for this purpose, in [3, 4] — the method 
of discrete elements. In [5], ice was modeled by an elastic plate lying on the surface of a stratified fluid. Models that 
allow cracks were considered in [6, 7]. Models of ice strengthened with reinforcing elements were presented in [8, 9]. 
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At this, in some papers, an ice cover is considered as an elastic plate lying on the surface of a reservoir [10, 11]. At 
the same time, in [12], on the basis of the conducted research, it is concluded that in some cases, the properties of ice are 
best described by the Kelvin-Voigt rheological model with one parameter (damping time). Therefore, numerous 
researchers use a viscoelastic plate when modeling the ice cover [13]. In [14], nonlinear models were used to describe the 
properties of ice. 

In some articles, the effect of a moving load on the ice cover was considered. In [15], the impact of a mobile load on 
the ice cover in a frozen channel was studied. In [16], the action of a load with an impulsive movement on the ice cover 
was described. Paper [17] is devoted to the study of the load moving along a frozen riverbed. Here, the rectilinear motion 
of the load was investigated [18]. However, in real conditions, it is often needed to deal with a load moving in a more 
complex way. Therefore, the objective of this work is to develop a method for solving problems about the action of a load 
moving over an ice cover along a complex trajectory. This will provide for a more accurate investigation of the effect of 
vehicles moving in a complex way on ice. 

This work is a continuation of research related to the problems about the effect of a moving load on various objects, 
whose results are presented in papers [19, 20]. 

Materials and Methods. Setting the task. A reservoir of finite depth with an infinite ice cover (an infinite plate), 
which is subject to the action of a vertical force moving in an arbitrary way — impulsively, is considered. It is assumed 
that the reservoir fluid is incompressible and executes irrotational motion. 

The problem is reduced to a system of differential equations [15]: 


x,y,t 
(14190, AW +07 O;W +kW +b0,F|_ = Olnrt) 
AF =0, 
where W(x,y,t) — ice cover deflection; E and 4 — Young's modulus and Poisson's ratio of ice, respectively; 


D = Eh3/12(1-7) — cylindrical bending stiffness; 4 — ice cover thickness; t. — strain relaxation time; A,” = (0,7+ 0,7)’; 
A= 0+ 0,?+0/; pxand ps— density of ice and water, respectively; c? = p,h/D; k = peg/D; b= ps/D; O(xy,t) — load 
acting on the ice surface; F(x,y,z,t) — speed potential. 

Under boundary conditions at z = 0 (ice-water boundary): 

OW =0.F. 
At the bottom of the reservoir at z = —H: 
0,F =0. 

In addition, it was assumed that the ice cover and the fluid in the reservoir were resting at the start time. The load was 
a concentrated unit force (one Newton) OQ (x, y, t), which moved arbitrarily along open free-form curve y. It was assumed 
that the displacement of the force was given in the form Q = Q(s(t)), where s — arc coordinate measured from some fixed 


= x(t) 


: x : 
point of trajectory y. The trajectory of movement was set parametrically in the form | ‘ where ¢ — time. 
= Volt 


The moving load was approximated by the expression: 
O(x,y,t) = 6° exp(-e" ((x —Xp (t))” +(y — Vo (1))” )\/e 


where  — numeric parameter. 
After applying the integral Fourier transform with respect to variables x and y, the integral Laplace transform with 
respect to ¢, formulas for calculating unknown functions W and F were obtained: 


W (x, y,t) mall 2p? /4s” Jo(pR(t-1))= ‘le ae O° — et" |dpdr, 


Pesan) asin [per tsa pete 


00 


to 
W (x,y,t)= xo) fren (pR(t -t)) 
00 


a [ete olata)s Japan, 
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Using the known relations from the theory of thin plates and the theory of the potential ideal flow, it is possible to 
obtain relations for calculating displacements and stresses in the ice cover, as well as the velocity vector component of 


fluid particles. 
When calculating the improper integral through numerical methods, the approximate relationship 


[, fle) did not 


io) A 
| f ( p) dp = | f ( p) dp was used, in which value A was chosen so large, that the error estimate 
0 0 


exceed the set value. 
Thus, for the amount of ice deflection 
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this estimate has the form: 
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Similar estimates can be obtained for other calculated quantities. These estimates were used to determine value A. 
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ef /40? 
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In the calculations carried out, value A was chosen such, that the estimate 


ten 

ssl Uo(t, p)dpdt| 0.001. 

When calculating the repeated integral, Simpson's quadrature formula (for variable t) and Chebyshev quadrature 
formula with equal weights for two nodes (for variable p) were used. The other values were calculated in the same way. 

Research Results. A method has been developed for solving problems on the action of a load moving along the ice 
cover of a reservoir filled with ideal fluid along a complex trajectory with variable speed. Using this method, calculations 
were carried out that showed the degree of impact of various parameters on the deformation of the ice cover. 

The described method does not impose restrictions on the shape of the trajectory of the concentrated force. In the 
calculations, a special case of a trajectory consisting of arcs of circles was considered (Fig. 1). The red dot indicates the 
position of the concentrated force at the time under consideration and the direction of movement of the force. 
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Fig. 1. Trajectory of concentrated force 
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The following parameter values were taken into account in the calculations: ice cover thickness h = 0.25 m, Young's 
modulus E = 500,000,000 H/m, Poisson's ratio of ice 1 = 1/3, ice density p = 900 kg/m’, fluid density p = 1,000 kg/m’, 
¢ =5. The calculation results are presented below. 

Figure 2 shows the change in the deflection of the ice cover at speed of force v = 2.5 m/s, tangential acceleration 


w; = | m/s”, reservoir depth H = 25 m, and relaxation time t) = 1s. 
The law of force motion along the trajectory was taken as: 
§=a,t? + aot? + ay. 
Coefficients a1, a2, a3 were selected in such a way that the force, being at the same point of the trajectory, had the 


required speed and tangential acceleration. 
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Fig. 2. Change in ice cover deflection 


At other values of these parameters, the qualitative nature of the distribution of ice cover deflections remained 
almost unchanged. 
Z,m 


Y,m 


Fig. 3. Fluid motion caused by action of moving load on the ice cover 
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The fluid motion caused by the action of a moving load at the same values of velocity, tangential acceleration of the 
load movement, relaxation time, and depth of the reservoir is shown in Figure 3 (the distribution of the velocity vector of 


fluid particles is shown). 
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b) 
Fig. 4. Change in amount of maximum deflection of ice cover depending on: 


a — tangential acceleration rate; b — depth of reservoir 


The effect of the tangential acceleration of the force movement on the maximum deflection of the ice cover is shown 
in Figure 4 a. In this case, the force speed was equal to v = 17.5 m/s, and the relaxation time 1) = | s. 

Figure 4 b shows a graph of the dependence of the maximum deflection of the ice cover W on the reservoir depth H. 
Here, the speed of the load movement was v = 17.5 m/s, tangential acceleration w; = 1 m/s”, and relaxation time T = Is. 
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Fig. 5. Change in maximum deflection of ice cover 
depending on force speed 


The dependence of the maximum deflection of the ice cover on the force speed is graphically shown in Figure 5. The 
depth of the reservoir was assumed to be 25 m, and the tangential acceleration — w,; = 1 m/s’. The solid line shows the 
dependence corresponding to the relaxation time To = | s, the dotted line corresponds to the relaxation time T = 10s. 
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Discussion and Conclusion. The influence of the reservoir depth on the maximum deflection of the ice was studied. 
A picture of the deflection of the ice cover was obtained due to the action of a concentrated force moving along a complex 
trajectory with variable speed. Calculations showed that with increasing depth of the reservoir, the maximum deflection 
of the ice cover decreased (Fig. 2). At the same time, a noticeable dependence of the deflection of the ice cover on the 
depth of the reservoir H occurred only for H<25 m. At great depths, the amount of the maximum deflections stabilized 
near a certain constant value and practically did not change. Thus, if H>25 m, then the depth of reservoirs can be 
considered infinite when calculating. 

An increase in tangential acceleration caused an increase in the deflection of the ice cover. Moreover, the dependence 
of deflection on tangential acceleration was very close to the linear dependence (Fig. 4). 

At low relaxation times tT, the speed of the load movement affected significantly the amount of ice deflection. At 
large times, the impact of the load movement speed on the deflection of the ice cover was noticeably reduced (Fig. 5). 

To study the state of the reservoir fluid, the distribution of the velocity vector of the fluid particles due to the action 
of the moving force on the ice was determined (Fig. 3). 

The developed method of solving problems and the results obtained with its help can be used in the construction of 
ice roads, design and construction of runways on ice. 
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Abstract 


Introduction. Computer modeling allows engineers to make valid design decisions by accurately assessing the thermal 
characteristics of design objects. The implementation of digital twin technology in the process of designing technical 
facilities is the current direction of scientific research and development. To do this, it is necessary to develop computer 
models whose accuracy meets the requirements for digital twins. However, the scientific literature does not widely present 
the results of research aimed at implementing digital twin technology in the design process. The general issues related to 
the use of digital twins in various industries are mainly considered. Therefore, the objective of this study was the 
development of a digital model and a comparative analysis of the accuracy of calculations of thermal characteristics of 
the design object. 

Materials and Methods. The main tool for conducting the research was the methodology proposed by the authors for 
developing a computer model of thermal characteristics for the implementation of digital twin technology. The numerical 
solution was implemented through constructing a thermal model for calculating the temperature field based on the finite 
element method in the ANSYS engineering analysis system from ANSYS, Inc. (USA). For the analytical solution, a 
computer model of thermal characteristics developed on the basis of the state-space method, implemented in the ANSYS 
Twin Builder module, was used. The state-space model was matched to the behavior of the original thermal model through 
approximating the transfer function to the stepwise response of the thermal load using the time domain vector 
approximation method. Verification of the constructed analytical model was carried out in the engineering calculation 
system MATLAB from the MathWorks company (USA). The research was carried out for a 400V machine model 
manufactured by NPO “Stankostroenie” LLC, Sterlitamak (Russia). 

Results. The developed digital model makes it possible to calculate the thermal characteristics of the design object with 
high accuracy. The results of the comparative analysis showed a high degree of correspondence between the values of 
thermal characteristics obtained using the proposed digital model and the results of numerical simulation. The maximum 
error in calculating thermal characteristics did not exceed 0.1°C. 

Discussion and Conclusion. Computer modeling that combines numerical calculation methods and a scientific approach based 
on digital twin technology, provides obtaining the result as close as possible to the results of experiments. The digital model 
proposed in the study is an effective solution, since it provides performing calculations to evaluate thermal characteristics in real 
time, which is one of the most important requirements for the implementation of digital twin technology. 


Keywords: digital twin, complex technical object, computer modeling, computer-aided design, temperature field, 
thermal characteristics 
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AHHOTalna 


Beedenue. KomupiotepHoe MoyesMpoBaHue NO3BOJIAeT HHDKeHepaM MIPHHUMaTb OOOCHOBAHHBIe MpOCKTHBIe peuleHHA 3a 
cueT TOUHOM OCHKM TeMJIOBbIX XapaKTepHCTHK OOBeEKTOB MpoeKTHpoBaHHA. AKTyaJIbHbIM HallpaBJIeHHeM Hay4HBbIx 
MCCeqOBAaHHH HM pa3spaOoTOK ABJIAeTCA peaH3alHA TEXHOJIOrHH UM(PpoOBbIX ABOMHUKOB B IIpolecce MmpoeKTHpoBaHHA 
TeXHHYeCKHX OOBeKTOB. JIM 3TOrTO HeOOXOAHMO pa3pa0aTbIBATb KOMIIbIOTEPHbIC MOM, TOUHOCTb KOTOPBIX 
COOTBETCTBYeT TpeOOBaHHAM, IIPeXbABJIACMBIM K UM(PpoBbIM ABOMHUKaM. OHako B Hay4HoH JMTepaType HeocTaTOUHO 
IIMpPOKO MpeCTaBJICHbI Pe3yJIbTAaTbI HCCC NOBaHHH, HalipaBJICHHbIXx Ha peasiM3al{Hio TEXHOJOFH WHPpoBbIX JBOMHUKOB 
B Ipoljecce npoekTupoBaHua. B ocCHOBHOM paCCMaTpHBaIOTCA OOMIMe BOMPOCHI, CBA3AHHbIe C IPHMeCHeEHHeM WHpoBbIXx 
J{BOMHHMKOB B pa3JIM4HbIX OTPacsAx MpOMbULeHHOcTH. IlosToMy WesbiO WaHHOTO HCcIeqOBaHHA ABHJIaCb pa3pa0oTKa 
WMppoBwoli MOAeIM U CpaBHUTeIbHbIM aHasIM3 TOUHOCTH PpacueTOB TeIMIOBbIX XAPaKTePHCTHK OOBEKTA MIpOCKTHPOBAHHA. 
Mamepuansi u memoodvi. B xKauecTBe OCHOBHOrO MHCTpyMeHTa Id MpOBeyeHHA UCCeqOBAaHHA BBICTyMaeT 
TIpeqOKeHHad ABTOPaMU MeTOANKa pa3spa0oTKH KOMIIBIOTePHOM MOC TeEIIOBbIX XapaKTePHCTHK JIA peaM3alHu 
TeXHOJOrHH WHPpoBLix ABOMHUKOB. UuceHHoe pelieHHe peaM30BaHO IyTeM MOCTpOeCHHA TeIWIOBOM MOJeIM IA 
pacueTa TeMiepaTypHOro MOJIA Ha OCHOBe MeTOa KOHCUHBIX 9JIEMEHTOB B CHCTeMe MH2KeHEpHOroO aHasiu3a «Ansys» OT 
kommaHun «Ansys Inc» (CLUA). Aa anammTuueckoro pelieHHa IpHMeHsAeTCA pa3pabOoTaHHad Ha OCHOBe MeTOsIa 
IIpOCcTpaHcTBa COCTOAHHM KOMIIbIOTepHad MOJ[elIb TEMIOBbIX XapaKTePHCTHK, peasIM30BaHHad B Mogyse «Ansys Twin 
Builder». Moyen mpoctpaHcTBa COCTOAHHM MpHBOAMTCA B COOTBETCTBHE C TOBeCHHeEM HCXOAHOM TeIWIOBOH MOesH 
IYTeM IIpHONMKeHHA MepexaTOuHON (yHKUMH K MOWaroBOMy OTKIMKy TemIOBOM Harpy3KH C IIPHMeHeHHeM MeTosa 
BeKTOPHOH alilpOKCHMalHH BO BpeMeHHOH OOsacTH. Bepudukalua NOCTpOeHHOM aHasMTHA4eCKOM MOJeJIM BBITIOHATACh 
B CHCTeMe HH>KeHepHBIX pacyeTos «MATLAB» ot komnanuu «The MathWorks» (CHIA). UccneqoBpanna npoBoqusuch 
Ia cTaHka Moyen 400V mpou3BozcTBa upexupuatua OOO «HIIO «Crankocrpoenne» r. Crepmutamak (Poccns). 
Pexynemamoi ucciedoeanua. PazspadoTaHa udpoBad MojelIb, WO3BOJIAFOWIad C BbICOKOM TOUHOCTHIO BBITIOHMTb pacueT 
TeIVIOBbIX XapaKTepHCTHK OOBeKTa TIpoeKTHpoBaHHA. PesyIbTAaTbI CpaBHUTebHOTO aHasIH3a TIOKA3bIBalOT BbICOKYIO CTeICHb 
COOTBETCTBHA 3HAYCHHM TEIMVIOBbIX XapaKTePHCTHK, MOJIYYCHHBIX C TMOMOMIbIO TIpeWIOKCHHOM WAPpoBol Moses, pesyIbTaTaM 
YMCIeHHOrO MoyesMpoBanHa. MakcuMasibHad MOrpelwHOCTb pacueta TeMIOBbIX XapaKTepHCTHK He peBbriaer 0,1°C. 
Ooécyotcodenue u 3akiio“enue. KomibroTepHoe MOeMpoBaHue, COUeTAIOee YHCJICHHbIC METObI pacdeTa U HaydHblii 
TOAXO, OCHOBaHHbIM Ha TeXHOJIOTHH UMpoBbIX JBOMHHKOB, MO3BOJIAIOT MOYYHTb pe3syIbTaT MaKCHMaJIbHO 
TIpHOJMKeHHBIM K pe3yIbTaTaM 9KcTIepHMeHToB. IIpeqoweHHad B YCCIeqOBaHHM WMdpoBat MOJeIb ABILACTCA 
93(pdbeKTHBHEIM PelIeHHeM, MOCKOJbKY MO3BOIACT BBIMIOJIHHTb pacueThI JIA OLCHKH TeIMJIOBbIX XAPaKTEPHCTHK B PerKHMe 
peasIbHOrO BpeMeHH, 4TO ABIACTCA ONHHUM U3 BarxkHeiWUx TpedoBaHHit IpH peaM3aljMu TeXHONOrHM WH@poBBIx 
JIBOMHHKOB. 


Ksnrouespie = cJIOBAa: IM @poBolrt JIBOMHUK, CIOKHBIM TeXHH4eckHh o0bexT, KOMIIBFOTepHOe MOeCIMpOBaHHe, 
aBTOMaTH3HpOBaHHOoe NpOCKTHPOBaHHe, TEMUepaTypHoe Noe, TCMIOBbIC XapaKTCPpHCTHUKU 


Baarojxapnocrn. ABTODBI BbIparxKaroT OaroqapHOcTE peak U peveH3eHTaM 3a BHHUMAaTeCJIbBHOeC OTHOMICHHE K CTAaTbe 
HW yKa3aHHble 3aMeCuaHHA, KOTOPble NO3BOJIMJIM NHOBbICHTb Ce KadeCcTBO. 


@unancuposanne. McciezoBaHve BbIMOJIHeHO Ip (PUHAHCoBOH NoszWep»xKKe MunuctepcTBa HaykKH H BBICLUeTO OOpa30BaHHA B 
paMkax TIporpaMMBI CTpaTerH4ecKoro akayleMH4eckoro sMepcTBa «I IpHoputet-2030» (cormamrenue Ne 075—15-2023-151). 


Aaa waTuposanna. [o3sesankun B.B., Nonaxos A.H. Peanus3ayua uHppoBol Moses TemIOBbIX XapaKTepHCTHK Ha 
ocHoBe TemMepaTypHoro oa. Advanced Engineering Research (Rostov-on-Don). 2024;24(2):178-189. 
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Introduction. Computer modeling has traditionally been an effective tool for solving thermal problems at an early 
stage of designing complex technical facilities. However, solving the problem requires an accurate assessment of the 
thermal characteristics of the design object to reduce the negative effects caused by an increase in temperature [1]. At the 
same time, one of the effective tools for preliminary assessment of thermal characteristics is simulation modeling in an 
engineering analysis system based on advanced digital solutions and developments. In [2], for example, a developed 
digital twin for determining thermal characteristics was presented by the authors Jianying Xiao and Kangoo Fan. The 
principle of operation of the twin was to simulate the thermal behavior of an object through displaying and correcting 
thermal boundary conditions. The experimental results showed a high accuracy of the model (more than 95%), which was 
essential for improving the accuracy of modeling thermal characteristics and thermal optimization. Therefore, the urgent 
direction of scientific research and development in the field of modeling is the use of artificial intelligence systems [3] 
and digital twins [4]. In [3] by Haoran Yi and others, an interactive model for correcting thermal boundary conditions 
based on a neural network was proposed to improve the accuracy of the analysis of thermal characteristics. The 
experimental results showed that the accuracy of calculating the temperature field exceeded 98%, and the accuracy of 
predicting thermal deformation was 96%, which effectively increased the simulation accuracy. In addition, in [5] by 
Kurganova N.V. and others, it was noted that digital twins were often used to improve physical prototypes of complex 
technical objects, since they not only allow for the information support for the design process, but also contribute to 
effective design decisions based on developments in the field of artificial intelligence. 

One of the characteristic features of digital twin technology is that reduced-order computer models are often used for 
simulation [6]. Therefore, the development of computer models is one of the basic conditions for the implementation of 
digital twin technology [7]. Model order reduction is an effective and mathematically understandable approach to 
overcome the time constraints of multidimensional simulation models. In [8] by Mirzoev D.A. and others, for example, 
a simple analytical model of thermal fields was proposed for the development of digital counterparts of the industrial arc 
welding process. Bordachev E.V. and Lapshin V.P. [9] presented the results of mathematical modeling of the temperature 
in the tool-product contact zone under metal turning. This approach provides obtaining an accurate assessment of thermal 
characteristics corresponding to the results of numerical experiments in real time. Schréder C. and Matthias V. [10] 
presented a reduced-order model and proposed a new model balancing procedure based on the transformation of the state 
shift. In conclusion, they presented the results of a comparative analysis and the results from literature sources obtained 
through a series of numerical experiments. 

The use of a linear and time-invariant reduced-order computer model allows for fast simulation while maintaining 
high calculation accuracy [11]. When developing a computer model, approximation [12] of the transfer function is 
performed to approximate the state space model to the step response of the initial thermal model [13]. Since the step 
response of the thermal load is derived from the base thermal model, the digital model should provide the same values of 
thermal characteristics. 

However, despite the fact that recently there has been an increase in interest in the digital twin, the scientific literature 
does not widely present the results of research related to the implementation of digital solutions in the design of technical 
facilities. Based on a systematic review of the literature and thematic analysis of publications on digital twins, one of the 
key knowledge gaps associated with the development of mathematical, software and methodological support for high- 
precision computer models in the framework of the implementation of digital twin technology has been identified. 

In this regard, the objective of the study was to develop a reduced-order computer model and analyze its feasibility as 
part of a digital model for accurate calculation of thermal characteristics of complex technical design objects. 

To achieve that, it was necessary to build a thermal model and calculate the temperature field of the design object, 
generate independent step responses of the thermal load using the developed software scenario [14], implement a digital 
model for calculating thermal characteristics, determine the error of calculations obtained using numerical and analytical 
solutions, and conduct a comparative analysis of the simulation results. 

Materials and Methods. The construction of the temperature field of the modeling object was performed for a 
homogeneous isotropic body on the basis of the equation of nonstationary thermal conductivity: 


6T/0t = a° AT +q,/c, (1) 
where 7 — temperature (°C); ¢— time (s); a= Jd/(p-e) — diffusivity coefficient (m7/s); 2 — thermal conductivity 
coefficient (W/m:°C); c — specific heat capacity (J/kg-°C); p — material density (kg/m*); A — Laplace operator; 
qv — volumetric heat dissipation power (W/m’). 
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The heat flow in the heat transfer process was assumed to be equal to the amount of heat transferred through an 
arbitrary surface area S per unit time t. It is expressed by the following equation: 


== J J ana. (2) 


where Q — heat flow (W); g, — heat-flux rate (W/m7); S — surface area (m7). 

The density of the heat flux during heat transfer was determined from the formula: 

qe =0°S-(Ts-Ta), (3) 

where a — heat transfer coefficient (W/(m°C)); Ts — surface temperature (°C); T.. — ambient temperature (°C). 

In this regard, to calculate the temperature field of the modeling object according to formula (1), heat flows (2) and (3) 
were set, which determined the amount of heat passing through the surface per unit time. 

A component (Fig. 1) of a metal-cutting machine 400V model manufactured by NPO “Stankostroenie” LLC 
(Sterlitamak, Russia) in the form of a drillhead was selected as the object of modeling. 


Fig. 1. Geometric model of the object in ANSYS Design Modeler: 1 — housing; 2 — electric motor; 
3 — spindle assembly; Q — heat flows; qv — volume heat release power; qn — heat flux densities; 
a — heat transfer coefficient; T1-T4 — temperature sensors 


Since the amount of heat released mainly by the electric motor (electromagnetic losses) and the spindle assembly 
(mechanical friction losses) was taken into account for calculating the temperature field, electric motor q,, and 
spindle assembly q,2, were taken as internal heat sources, near which the corresponding heat flows Q; and Q2 were 
set. The densities of heat fluxes qn, qn2 were assigned to surfaces located near the electric motor and spindle 
assembly, respectively. 

Convection was determined by the heat transfer coefficient a taking into account the conditions of heat transfer 
(natural convection in air). Since the machine is in contact with a gaseous medium (air), the amount of heat given by the 
heated surface to the environment per time unit f, is directly proportional to the difference in temperature between surface 
Ts and medium T., depending on the area of the heat-emitting surface S (3). 

When constructing a thermal model of an object (solid) consisting of a homogeneous material (structural steel) with 
constant thermophysical properties and the presence of internal heat sources, the following initial and boundary conditions 
were assigned: 

— initial conditions took into account the fixation of a constant temperature over the entire surface of the modeling 
object (¢ = 0: T = To = const); 

— boundary conditions of the second kind were set by the heat flows of the electric motor (Q)), the spindle assembly (Q2), 
the density of the heat flux (gn1) from the electric motor to the front wall of the housing and (q,2) to the inner surfaces; 

— boundary conditions of the third kind were set by the heat transfer coefficient (a) for the surfaces located inside the 
body of the drillhead; 

— boundary conditions of the fourth kind for the contact joints of the surfaces took into account the perfect thermal 
contact and the absence of thermal resistance: 

oT 
"én 

Heat flows (Q:, Q2), heat flux densities (gn, qnz), aS well as volume heat dissipation capacity (qui, gz), Were assigned 
in accordance with well-known recommendations for metal-cutting machines [1]. Boundary conditions, as well as heat 
flows, were set for external surfaces. Therefore, it was taken into account that the interrelationships of thermal fields were 


T| 5 =T aid +g: (4) 


+0’ 


a on 


present only between the outer surfaces. 
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Differential equation (1), together with the initial and boundary conditions of the second (2), third (3) and fourth (4) kinds, 
is a mathematical formulation of the problem. The task was solved using numerical and analytical modeling methods. 

The numerical solution was performed on the basis of the finite element method in the ANSYS engineering analysis 
system, which is being developed by ANSYS Inc. (USA) and supplied by “Modeling and Digital Twins” AO, an 
authorized ANSYS distributor in Russia. The geometric model (Fig. 1) of the object was imported into the ANSYS 
Workbench project with the addition of the Transient Thermal analysis block. In the numerical solution, the simulation 
parameters and boundary conditions of the thermal model were calibrated (Table 1) to approximate the model temperature 
values to the experimental data. 

Table | 
Boundary conditions (parameters) of the thermal model 


issipati i “as Heat t fe 
Heat dissipation capacity Heat flows Heat flux densities ee 


Parameter coefficient 
qui, Wim? qu2, Wim? O1, W O2, W gu, W/m? gu2, W/m? a, W/(m?-°C) 
Value 6,500 1,000 28 15 32 18 15 


In the ANSYS Mechanical module, the initial (initial temperature Tp = 24 °C) and boundary (Table 1) conditions were 
assigned for the developed grid model, the temperature field was constructed (Fig. 2). In this case, the thermal 
conductivity coefficient was assumed to be equal to 4 = 60.5 W/(m-:°C) and was assigned as such for structural steel. 

The thermal model of the object included two contact connections for an electric motor and a spindle cartridge with a 
drillhead body, and contained 7 thermal boundary conditions. The developed grid model consisted of 16,309 elements 
and 58,527 nodes. 

The total simulation time of 21,600 seconds (6 hours) was divided into intervals (1 hour) and steps At = 360 seconds 
(6 minutes), a total of N = 60 steps within which the parameters of the thermal model (boundary conditions) were assumed 
to be constant and independent of time. 


Model (84) 
(4, Geometry 
83 Materials 
vik Coordinate Systems 
# te) Connections 
f--/) Mesh 
=) [4 Transient Thermal (85) 
v=o Initial Temperature 
lil] Analysis Settings 
2} Simplorer 
v®, Internal Heat Generation1 
y@, Internal Heat Generation2 
v&, Heat Flow1 
Vv, Heat Flow2 
JB Heat Flux1 
JB4 Heat Flux2 
vy Convection 
= Solution (B6) 
--/ 3} Solution Information 
y@® Temperature 
=} Simplorer 
ye TempProbe 1 
y® TempProbe2 
y%@ TempProbe3 
ye TempProbe4 


a 


Fig. 2. Temperature field of the object in ANSYS Mechanical 


The calculation of the temperature field was required to generate an independent step response of the thermal load, 
containing information about temperature variation over time (thermal characteristics) for each parameter (Table 1) of 
the thermal model separately: 

y:(t)=7,,(¢),@=Lk fel, (5) 
where A — number of temperature sensors; m — number of parameters of the thermal model. 

Step response generation was performed in the ANSYS Mechanical module using the Application Customization 
Toolkit (ACT) extension, which supported the implementation of user scenarios. This made it possible to develop a 
software script in the Python programming language to automatically generate a special set of files containing independent 
step responses [14]. 
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The analytical solution was performed on the basis of the state-space method through constructing a model of thermal 
characteristics according to the formula: 
x = Ax+ Bu 
. (6) 
y=Cx+Du 
where x = Ox / Ot — derivative of state vector x over time ¢; y and u — vectors of output and input data, respectively; 
A, B, C, D— matrices of constant coefficients. 


In equation (6), vector x=(x1,%,...,.xy)" contains state variables, input —_- vector 
U = (Gris Gv25 Q1, Q2, Anis Yn2, To1, To2, ...5 Tox)’ — values of the parameters of the thermal model (boundary and initial 
conditions), output vector y =(T\, T>, ..., Tx)’ — values of thermal characteristics. 


The transfer function model is expressed by the following equation: 


y,(t)= ES, (t)u,,i=Lk 
where H — matrix complex transfer function. 
Matrix transfer function H is obtained through applying the Laplace transform to formula (6), which is expressed by 


(7) 


the following equation: 
H(s)=C(sI—A)' B+D, (8) 
where s — Laplace complex variable; / — unity diagonal matrix. 

Transfer function (8) reflects the dependence of the Laplace transform of output variable Y(s) = H(s)U(s) on the 
Laplace transform of the input variable U(s) of model (6) under zero initial conditions x(0) = xo = 0. In this case, the 
dimension of the transfer function matrix H depends on the output value & = 4 and the dimension of the input value 
m = 10, which corresponds to the dimension of the initial step response (5). Therefore, model (6) is brought into line with 
the behavior of the source thermal model through approximating its transfer function (7) to step response (5) using the 
vector approximation method [15]. 

To construct the coefficient matrices of model (6), a reverse transition is performed from transfer function (8) to the 
model in the state space. In this case, transfer function (8) takes the form of the equation, whose denominator contains a 
characteristic polynomial of degree / = 4 (the order of the system), and the numerator contains a polynomial of degree 
z=1-1: 

ra 


G(s)= a. P(s)=bo + > as O(s)= aq + >, 28 (9) 


where a and b— coefficients of polynomials Q(s) and P(s), respectively. 

The roots of polynomials Q(s) and P(s) represent the poles and zeros of transfer function (8), respectively. The method 
of indefinite multipliers is applied to equation (9) to decompose each element of matrix H into elementary fractions. 
Denoting the poles of the characteristic polynomial by p;, we obtain an equation of the following form: 


G(s)=D+ 0" A (10) 


? 
i=l 8 — p; 


P(s) 


where R; = lim (s— p;)——~ — matrix of dimension (k x m); q — number of poles. 
Pi O(s) 
The rank of matrix R; is denoted by 7;, and its decomposition into the product of two matrices with the full rank of the 
column and row, respectively, is performed: 


REC BP” rank (R=, (11) 
The matrices of model (6) are diagonal, of dimension 4””, B’’”, C®", D®™ (n = 56) and contain elements that are 


obtained directly from the coefficients of transfer function (10). 
System matrix A and control matrix B contain the following coefficients: 


fp, 0 0 .. 07 b, O .. 0 
we bis 0 
0 0 0 b 0 
A= Pi2 : B- ql (12) 
0 0 PB 0 0 by Dim 
[0 0 0 Pm| [0 0 be 
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Output matrix C and feed forward matrix D contain the following coefficients: 


l ci wie Cig 0 eee 0 0 ane 0 0. ... 1 .. O 0 
x eee eee eee > = ose ,D= eee wee eae eae eae eae : 13 
0... O C(k-1)i C(k-1)q O.. 0 Oo .. 0 .. 1 0 ( ) 
| 0 eis 0 0 ees 0 Chie Chg 0... 0 a. 0 1 


This approach makes it possible to perform the transition from a model in the form of transfer function (8) to a model 
in state space (6). The transition algorithms are also implemented in the MATLAB engineering calculation system of The 
MathWorks, Inc. (USA), in the form of special functions «ss2tf()» and «tf2ss()». 

To obtain the values of thermal characteristics, the Cauchy problem is solved for a system of ordinary differential 
equations, since in formula (6), variable x is a derivative of the vector of states of the temperature field in time ¢. The 
solution to system of equations (6) is obtained using the fourth-order Runge-Kutta method. 

Verification of the constructed analytical model was performed through conducting a series of computational 
experiments using an application program developed in the MATLAB system, which included the implementation of the 
fourth-order Runge-Kutta method (Fig. 3). Computational experiments were conducted on a personal computer (AMD 
Ryzen 5 5600U processor with Radeon Graphics 2.30 GHz, RAM 16.0 GB, system type 64-bit Windows 10 Pro version 
21H2 operating system), whose characteristics were basic for modern computing technology. 

The constructed analytical model and the application of the fourth-order Runge-Kutta method for solving a 
system of differential equations made it possible to calculate the values of thermal characteristics with high accuracy 
(Fig. 3). Due to the fact that the maximum error values, i.e., the difference between the model values of thermal 
characteristics obtained using numerical (FE-Model) and analytical (LTI-Model) solutions, did not exceed 0.72 °C 
over the entire modeling interval. 
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Fig. 3. Graphs of thermal characteristics in MATLAB system 


Linear and Time-Invariant (LTT) Reduced Order Model (ROM) was developed in the ANSYS Twin Builder module 
of the ANSYS engineering analysis system based on the object's temperature field and the step response generated in the 
ANSYS Mechanical module. 

The implementation of the digital model of thermal characteristics based on a temperature field using vector 
approximation consists in the sequential execution of seven main stages. 
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Stage 1. Importing a geometric model of an object and building a thermal model in an engineering analysis system. 

Stage 2. Calculation of the object's temperature field based on the developed thermal model. 

Stage 3. Generation of an independent step response based on the results of numerical modeling of the thermal 
characteristics of the object. 

Step 4. Application of the vector approximation algorithm to obtain the poles and zeros of the transfer function of the 
state space model. 

Stage 5. Building a model of the state space based on a known transfer function. 

Stage 6. Development of a computer model of thermal characteristics. 

Stage 7. Implementation of the digital model of the object. 

Thus, the proposed digital model containing a computer model for an accurate assessment of the thermal 
characteristics of complex technical design objects was obtained through sequential completing all the above steps. 

Research Results. The developed computer model (Thermal_SG400V_SML1) contains 6 inputs and 4 outputs (Fig. 4); 
this provides identifying the relationship between the parameters of the thermal model and the values of thermal 
characteristics. The volume heat release power (qv, g.2), heat flows (Q:, Qs) and heat flux densities (gn, qn2) were taken as 
input data of the computer model. The output data were thermal characteristics (71-74). The computer model is part of the 
digital model (Fig. 4 a) of thermal characteristics implemented in the ANSYS Twin Builder module. 

The values of the parameters of the digital model, presented in the tabular form (Fig. 4 c), are fed to the input of the 
computer model (Fig. 4 a) using the STEP components. The graphic module shows the values of thermal characteristics 
(Fig. 4 5) obtained at the output of the computer model. 


Thermal SG400V_ SML 


1°C 

34. 

34. 

29.0 

24. 

0 5,000 10,000 ~—-15,000 t, 8 
b) 
1 2 3 4 5 6 
Time, s 3,600 7,200 10,800 14,400 18,000 21,600 

qn2 18.0 18.0 18.0 18.0 18.0 18.0 
qvl 6,500.0 6,500.0 6,500.0 6,500.0 6,500.0 6,500.0 
Ql 28.0 28.0 28.0 28.0 28.0 28.0 
Q2 15.0 15.0 15.0 15.0 15.0 15.0 
Ov2 1,000.0 1,000.0 1,000.0 1,000.0 1,000.0 1,000.0 
qn\ 32.0 32.0 32.0 32.0 32.0 32.0 


Fig. 4. Implementation of digital model in Ansys Twin Builder system: a — computer model of thermal characteristics; 
b — module for graphical representation of results; c — module for monitoring input data of the computer model 


During the development of the computer model, the limits of dimension from 2 to 4 orders, values of the target error 
¢ = 5x10°, and the tolerance for the zero order ¢) = 2x10 were set. The remaining parameters were set automatically, 
since the vector approximation method was as automated as possible in comparison to other methods that were supported 
in the ANSYS Twin Builder module. 

In the process of developing a computer model, the ANSYS Twin Builder module automatically generated a special 


matrix of approximation errors M ; =|| y;(t)- y ;(¢)||; in the time domain, each element of which reflected the difference 


between the values of the thermal characteristics of step response (5) and transfer function (7) of the model in the state 
space. It is expressed by the following equation: 
0.43 4.04 1.53 4.97 2.22 1.82 
0.21 2.22 2.54 1.07 1.16 1.95 
M= x10, (14) 
2.66 1.63 0.39 0.26 3.43 1.82 


1.46 3.21 0.06 0.84 2.95 3.50 
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In this case, the maximum relative error did not exceed value ¢ = 4.97 x 10°. At the same time, all other error values 
turned out to be less than the specified limit ¢ = 5 x 10>. Zero-error value in matrix (14) meant that the input was ignored 
due to a very small contribution. 

The estimation of the accuracy of the calculation of thermal characteristics using the proposed digital model was 
carried out through comparative analysis of the results obtained using numerical and analytical solutions. The comparative 
analysis was performed according to the criterion of maximum error. Maximum error ATinax, i.¢., the difference between 
the values of thermal characteristics obtained using numerical and analytical solutions for all temperature sensors, was 
calculated at each time by the formula: 

AT, 


max 


= max AT, (15) 


jalan 
where AT; = |Zj— Tj,<| — error; T,; and T;,— temperature values of the finite element and digital models, respectively 
G= 1, m); m — number of temperature sensors. 


To assess the digital model accuracy, an error calculation was performed, whose results were presented in the form of 
a surface (Fig. 5 a) and a linear graph (Fig. 5 5) of the maximum time error. The surface (Fig. 5 a) represented the 
calculated error values for each temperature sensor (dT1, dT2, dT3, dT4) individually and at each time point over the 
entire simulation interval. 
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Fig. 5. Errors in the calculation of thermal characteristics: a — error for each temperature sensor; 
b — maximum error for all temperature sensors 
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The graph (Fig. 5 b) showed that for the selected time interval, maximum errorA7ina, did not exceed 0.1 °C. The error 
in calculating the thermal characteristics for each temperature sensor separately (Fig. 5 a) did not exceed the specified 
limits either. 

Discussion and Conclusion. The obtained results of computational experiments and the conducted comparative 
analysis confirm the efficiency of the proposed digital model, which provides calculating thermal characteristics with 
high accuracy ATina. = 0.052 °C, for further analysis and identification of the thermal model. 

The temperature field of the design object is influenced by many factors, which complicates the determination 
of the nominal values of thermal boundary conditions. To solve this problem, the study first analyzes the key 
technologies involved in the implementation of the digital twin of a complex technical object, then builds a thermal 
model and a computer reduced-order model LTI ROM of the transient temperature field. The developed computer 
model is used as part of a digital model that provides obtaining an accurate assessment of the thermal characteristics 
of the design object and thereby increases the efficiency of design procedures in the process of computer-aided 
design of complex technical facilities. 

The accuracy and efficiency of calculations of the reduced-order computer model and the source full-order thermal 
model are evaluated by comparative analysis of the simulation results. The results of computational experiments show 
that, from the point of view of calculation accuracy, computer models of reduced order and finite element models of full 
order are generally comparable in accuracy, the maximum calculation error is within the acceptable range and does not 
exceed 0.1°C. 

The results obtained during computational experiments do not contradict the results presented in the sources of 
scientific literature on similar topics and allow us to conclude that the use of the proposed digital model is effective for 
evaluating the thermal characteristics of complex technical objects in real time, which is one of the most important 
conditions for the implementation of digital twin technology. 

However, changes in the temperature field of a complex technical facility still depend on many factors. Therefore, in 
further research, it is necessary to develop computer models of thermal deformations and conduct effective optimization 


algorithms based on artificial intelligence to provide the reliability of simulation results obtained using digital twins. 
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Abstract 


Introduction. The dynamic loads during the start-up of a bridge crane can cause excessive stress in the structure and 
components, leading to potential safety hazards and increased wear and tear. To reduce the influence of the dynamic 
loads, various strategies can be implemented including optimization of the acceleration and deceleration profiles, using 
the soft start controls, implementing the vibration damping systems. It is vital to ensure that the proper crane maintenance 
and inspection protocols are in place. By reducing the impact of dynamic loads during the start-up, the overall 
performance and longevity of a bridge crane can be improved, ultimately enhancing safety and efficiency of the industrial 
operations. The present research offers a new approach to improving the efficiency and safety of industrial operations by 
providing a more precise account of the dynamic loads during the start-up of a bridge crane. The objective of this study 
is to develop a mathematical model for investigating the mechanical properties of the bridge cranes by analyzing the 
dynamic loads that occur during lifting operations. 

Materials and Methods. The development of the mathematical model was based on the kinetic model of the system, 
which included three connecting blocks and two flexible connections for a more accurate description of the bridge crane 
structure. Lagrange’s equations incorporating the information about the geometry and structure of a bridge crane were 
used. They made it possible to describe the motion of a system with the multiple elements and degrees of freedom. 
Processing and analysis of the results of the mathematical model were carried out in the MATLAB program using the 
Runge-Kutta method. 

Results. As a result of the research, a mathematical model was developed to study the dynamic loads affecting a bridge 
crane during lifting operations. Graphs describing the dependences of speed, acceleration, load, and rope angle over time, 
and their influence on the crane beam were plotted. The changes in these parameters over time, including their maximum 
values, were analyzed. The reasons for load changes and factors influencing the extension of lifting machines’ service 
life as well as reducing metal consumption during production thereof were identified. 

Discussion and Conclusion. The developed mathematical model and its numerical solution using the specialized software 
(MATLAB) allow for conducting the dynamic analysis of the bridge crane structures and determining the optimal design 
solutions. The analysis of the factors influencing the load changes leads to the conclusion that the use of this model can 
significantly reduce the load magnitudes and metal consumption, as well as increase the service life of lifting machines. 
The results obtained with the developed mathematical model and its numerical solution are useful for optimizing the crane 
structures, providing compliance with the operational requirements, and extending the service life of lifting machines. 


Keywords: bridge crane, dynamic load, kinetic model, load lifting, dynamic analysis 
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AHHOTalna 

Beedenue. unamuueckne Harpy3KH BO BpeMa 3allyCKa MOCTOBOTO KpaHa MOTYT BBI3bIBaTb H30bITOUHBIe HallpsoKeHHA B 
KOHCTPYKUMH, IPWBOAA K MOTCHIMMAIbHBIM PHCKaM WM yBemM4eHuIO u3HOCA. Jd CHWOKeHHA BIIMAHHA JWHAMMUeCKHX 
Harpy30K MO%XKHO IIPHMCHATh pa3IMYHble CTpaTeruu, BKIHOUad ONTHUMUM3allMIo Wpodusleti ycKkopeHua U 3aMeyIeHua, 
MCHOb30BaHHe IIaBHOTO IycKa, BHeypeHue cucTeM amMopTu3alnNH. BaxxkHo oOecieduBaTb MCHOHeHHe paBMJIbHBIX 
TMpOTOKOJIOB OOcCIyKMBAHHA HW MHCIeKWMH KpaHos. IlyTem cHwKeHUA BO3TeHCTBUA JMHAMMYeCKHX Harpy30K BO BpeMa 
3allyCKa MO2KHO YJIVUINMTb OOIMLYIO MPpOH3BOAMTeIbHOCTh HW JOMTOBCYHOCTb MOCTOBOFO KpaHa, MOBbICHB B KOHCUHOM 
uTore O6e30nacHocTs u 3deKTHBHOCTh MPOMBIMINICHHBIX OlepalHi. JjanHoe uccieqOBaHue MpesaraeT HOBBIM NOAKoT 
K MOBBIINeCHHIO 3P@eKTHBHOCTH WM Oe30MaCHOCTH MPOMBINWICHHEIX OMepaluii 3a cueT Oomee TOUHOTO yueTa 
MHaMMYecKuX Harpy30K MOCTOBOTO KpaHa pu mycke. Lem padoTEr — pa3padoTKa MaTeMaTHYeCcKOM MOeIM TA 
W3y4eHHA MCXaHHYeCKHX CBOMCTB MOCTOBBIX KpaHOB ITyTeM aHasI3a JMHAaMMYeCKHX Harpy30K, BO3HHKAIOINHX BO BpeMa 
TIOJbCMHBIX Olepalluii. 

Mamepuanot u memooot. PaspadoTka MaTeMaTH4ecKOH Moe Oba BINOWHeHa Ha OCHOBe KHHeTH4eCKOH MOE 
CHCTeMBI, BKJIOUAIOIeH Tp COCIHHUTeIbHEIX OOKa HM Ba TMOKHX coeqMHeHUA WIA OoMee TOUHOTO onmucaHHa 
KOHCTpyKUMM MocToBoro KpaHa. Vcromb30BaHbI ypaBHeHua JlarpaHoxa, BKIHOUAONIMe HHDOpMalluIo O TeoMeTpHN U 
CTpyKType MOCTOBOrO KpaHa. OHM NO3BOIMIM OMMCaTb BIOKeHHe CUCTeMBI C HECKOJIBKUMH 39JICEMCHTaMH UM HECKOJIBKMMM 
cTeleHaMu cBoOorBI. OOpaboTka H aHasIM3 pe3yIbTAaTOB MaTeMaTHUYeCKOM MOJeIM ObLIM MpOM3Be{eHbI B IporpaMMe 
MATLAB c npuMeHenuem Metoya Pyure-KyTtsl. 

Pezyjemamol ucciedoeanua. B pe3ynbTaTe uccileqoBaHHua Oba pa3spaboTaHa MaTeMaTHYecKad MOJ[elIb JIA H3y4eHuA 
WHaMM4ecKHXx Harpy30K Ha MOCTOBOM KpaH BO BPeMA MOJbeMHBIX Oepalul. ToctpoeHsi rpaduku, onMcHrBarolue 
3aBMCHMOCTH CKOPOCTH, yCKOpeHHa, Harpy3KH WM yrla KaHaTa OTHOCHTeCIBHO BPCMeHH VM UX BIMAHMe Ha Oalky KpaHa. 
IIpoaHamM3upoBaHo W3MeHeHHe 3THX WapaMeTpoOB BO BPeMeHH, BKJIKOUaT HX MaKCHMAJIbHbIe 3HaYeHHA. OrlpeseuIeHbI 
IIPHUMHBI H3MCHCHHM Harpy3KH MW *aKkTOPHI, BIIMAFOMIHe Ha yBeIM4YeHHe CpoKa CIy2KObI HM CHWKCHME MeTAaIIOCMKOCTU 
IIPH MIpOW3BO]{CTBe TIOXbBeCMHBIX MallIHH. 

O@6cystcoenue u 3aknto“enue. Paspad0oTaHHad MaTeMaTW4eCcKad MOJeIb U Ce YHMCMCHHOe pellleHve C MCIOIb30BaHHeM 
cilel{HasiM3HpOBaHHOro MporpaMMHoro obecneyenHua (porpammMa MATLAB) no3Bous10T NpOBOAMTh AHHamMMyeckHit 
aHalIv3 KOHCTIpyKUMM MOCTOBOTO KpaHa HM OMpeeATb OMTHMAaIIbHbIC KOHCTPYKTHBHbIe pellienua. AHaiu3 dakTopos, 
BIIMAIOMWIMX Ha W3MCHeHHe Harpy3KH, MO3BOIAeT CieaTb BEIBOT, YTO IPH UCMOb30BaHHM WaHHOM MOJesIM MOKHO 
3HAYHMTeIbHO CHH3UTb BeJIMYHHY Harpy30K HM MeTAJIIOCMKOCTE, a TalOKe YBCJIMYUTb CpOK CIIy2KObI WOJbCMHBIX MallIHH. 
Pe3yibTaTbI, MOUYYeHHbIe Ip MOMOMM paspaboTaHHOM MaTeMaTH4eCcKOM MOJelIM, H ee UMCIICHHOe pellieHHe Mose3HBI 
IIPH OMTHMU3allM KOHCTpYKIMH KpaHosB, OOecieyveHHH COOTBETCTBHA OllepallHOHHBIX TpeOoBaHHi U MpomeHuN cpoKka 
CJLY2KObI WOUbEMHBIX MallIHH. 


KunroueBble CJI0Ba: MOCTOBOM KpaH, THHaMvYecKad Hal py3Ka, KHHCTHYICCKAA MOJICJIb, TOBEM Tpy30B, TMHaMuueckHit aHasIn3 


lia waTupopannsa. AutnOac U.P. Moyemmpopanue WuHaMMyecKux Harpy30K, BO3/[eHCTBYIOWJMX Ha MOCTOBOM KpaH B 
Moment ItycKa. Advanced Engineering Research (Rostov-on-Don). 2024;24(2):190-197. https://doi.org/10.23947/2687- 
1653-2024-24-2-190-197 


Introduction. The analysis of dynamic processes in the mechanical part plays an important role in the development 
of new overhead cranes and modernization of existing ones in order to reduce loads on control devices and extend their 
service life [1]. 

The bridge crane is subjected to dynamic loads during non-static operations, such as acceleration and braking. Analysis 
of these processes allows identifying hidden impacts on the dynamic behavior of the bridge crane. Therefore, it is 
paramount for the researcher to make an optimal design choice to reduce these loads, ensuring that the crane can meet the 
required operating conditions [2]. 

In [3], a dynamic model of a crane lifting system was developed, using which an accurate direct numerical integration 
method was proposed for calculating the dynamic loads of the system. 
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Papers [4, 5] studied dynamic loads in a metal structure, taking into account fatigue of the metal material. However, 
the researcher neglected the impact of forces from the drive of the lifting mechanism operating with artificial parameters. 

In [6], dynamic loads in a bridge crane were determined during the operation of a lifting mechanism when hoisting a 
load suspended on a rope. The most important case studied was the effect of dynamic loads on the crane when removing 
a load from a solid foundation, at the moment of lifting-off. 

Steel structures of bridge cranes experience non-stationary loads with different stress amplitudes and asymmetry of 
the working cycle [7]. To study the real load of bridge cranes under typical operating conditions, constant recording of 
their stress state is required, which is labor-intensive [8]. Therefore, statistical processing of the obtained results is used 
to assess the loading elements of metal structures of bridge cranes [9]. This involves changing individual components of 
the total load of metal elements, such as the gravity of the load being lifted, the angle of rotation of the load, and weather 
loads, and then summing them according to the laws of probability theory [9, 10]. This approach is less labor-intensive 
than a comprehensive study of loading under typical operating conditions. However, determining the probabilistic 
characteristics of individual random loads also takes time, so the method for calculating load combinations is widely used 
in crane construction [11]. 

In [12], it was found that during the stage of selecting the rope slack, the value of the stator current of the electric 
motor of the lifting mechanism did not depend on the mass of the suspended load. However, as the load increased, the 
time of its rise also increased, and at the stage of separating the load from the surface, the amplitude values of the current 
increased. Furthermore, a noticeable difference appeared after five periods of mains voltage from the beginning of the 
stage. However, the researcher neglected the significant influence of forces arising from the drive of the lifting mechanism 
working with artificial elements. 

By using the developed mathematical model, the research aims at optimizing the design of bridge cranes through 
studying the dynamic loads that occur during lifting operations, ensuring that the crane has the ability to meet the required 
operating conditions. 

Materials and Methods 

Kinetic Model of the System under Study. When developing a kinetic model of the system under study, it can be 
represented that the construction of a bridge crane for a given motion form consists of three connecting blocks and two 
flexible joints, and has the form shown in Figure 1: 


Fig. 1. Bridge crane kinetic model 


mg — mass of the main bridge, moving along the X-axis; mp — mass of the lifting mechanism, moving along the Y-axis; 
mg — payload; Kz and Cg — elasticity modulus and damping coefficient of the main bridge; Kr and Cr — elasticity 
modulus and damping coefficient of the ropes; o — swing angle of lifting mechanism winding. 

Derivation of Kinetic Equations Reflecting the Motion of the Dynamic Model 

Mathematical equations representing the motion of the dynamic model are derived from the partial differential 
Lagrange equation, which is considered to be one of the best methods used specifically in cases where the system consists 
of more than one element and when it has several degrees of freedom. 
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where 7 — kinetic energy; U — potential energy; D — damping energy; Q — external forces affecting the whole system; 
qi — common system of coordinates; i — degrees of freedom of the model under study. 
Kinetic energy equation: 


T =0.5m,22 +0.5mpz2 +0.5mgzz +0.57-a”. (4) 
Potential energy equation: 

U =0.5K 323 +0.5K,(zg—-ra+z,) . (5) 

Damping equation: 

2 

D =0.5C 523 +0.5C, (2% -ra+zZ,) . (6) 

System equations: 
(mg +mp)2g+Kg2g+K, (2, -ratzg)+Cg2g+Cr(z,-ratig)=0 (7) 
‘mote +K,(zo-ra+zg)+C, [tg -ra+zg )=0 (8) 
"J-G&+K,r(a-—Zzg —Zg)+C,-r(a—Zg -23)=Ma-M,, (9) 


where Mz, M.— resistance momentum and lifting mechanism momentum; J — inertia of rotating mass of the lifting mechanism. 

Numerical Solution of a Mathematical Model Using the Fourth-Degree Equation of Runge-Kutta 

The numerical solution to equations (7-9) was obtained using the Runge-Kutta method in the MATLAB program. 
The dynamic equations were derived within the program, incorporating the input data and a set of commands to process 
these equations. The resulting graphs illustrate the interconnections between the various blocks and components of the 
crane structure under consideration. 

The Conditioning of the Input Values Required to Solve the Model 

The study was conducted on an ACE type bridge crane consisting of three parts (Fig. 2): 

— lifting trolley; 

— main bridge, which supports the lifting mechanism; 

— end trucks, which support the main bridge. 


Fig. 2. Main parts of bridge crane ACE 


To align the operation of the bridge crane with a standard set of coordinate axes, the following assumption was made: 

— the lifting trolley functions as a unit responsible for raising a load along the Z-axis and allows for horizontal motion 
along the Y-axis relative to the main bridge, which is the axis along which the crane moves. 

A study was conducted on a prototype bridge crane with a lifting capacity of 10 tons and a width of 21.5 meters. The 
following characteristics were considered: 

— payload: mg = 10,000 kg; 

— mass of the two main bridges: mg = 8,100 kg; 

— mass of the lifting mechanism with the trolley: mp = 700 kg. 

When determining the input values, all the laws of designing the structures of lifting devices were followed, and the 
connections of all components were taken into account. The main factors considered were: 

— power of the drive of the lifting mechanism along the Z-axis; 

— stiffness coefficient of the steel structure (Ka); 
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— stiffness coefficient of the rope (K,); 
— damping coefficient of the metal frame (Ca); 
— damping coefficient of the ropes (C,). 
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Table | 
Shows the elements of the design structure of the bridge crane that was studied. 

Element name Value Unit of measurement Notation 
Load 10 Ton mG 
Crane mass 8,100 Kg mB 
Trolley mass 700 Kg mp 
Crane length 21.5 m L 
Height of lift 5 m A 
Beam device degree 2 - A 
Gear box ratio 4.5 - im 
Coil radius 0.25 m R 
Rope diameter 16.5 mm dk 
Engine power 30 kW Nn 
Speed of the engine rotor core 905 r.p.m Mn 
Rope stiffness coefficient 11,169.8 N/mm Kz, 
Rope damping coefficient 23,934.4 N/mm K, 
Rope damping coefficient 83.37 N.sec/m oF 
Metal frame damping coefficient 30.29 N.sec/m Cp 


Research Results 


Study of the Model Operation Using a Computer 
Structural Behavior of a Metal Bridge Crane under Momenta Loads 
The structural behavior of a metal bridge crane during the process of lifting a load is depicted by three curves (Fig. 3). 
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Fig. 3. Coordinates of vertical motion, speed and acceleration of the COG of the bridge cranes versus time curves 


In Figure 3, the first curve shows the coordinates of the bridge crane along the ordinate axis as a function of time. 
During lifting, some vibration of the crane's metal structure is observed for 5—10 seconds, which then stabilizes and does 


not affect its rigidity. 


The second curve in this figure depicts the change in the speed of the center of gravity of the bridge crane over time. 
During lifting the load, the speed initially increases and then gradually decreases until stabilization. This indicates that 
the center of gravity of the bridge crane assumes a stable position within the specified time period, during which the 


vibration stabilizes. 


The third curve reflects the change in the center of gravity of the bridge crane over time. It is noteworthy that the 
acceleration value at the moment of lifting the load is 0.07 m/s?, with the dynamic load reaching its maximum.° 
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Load Curves 
Figure 4 shows three curves that reflect motion of the load during the crane operation at the time of lifting. 


Coordinate zB, m 


Speed dzB, m/s 


Acceleration zB, m/s? 


Fig. 4. Coordinates of the speed and acceleration of the load lifting versus time curves 


In Figure 4, the first curve reflects the change in the position of the center of gravity of the load over time during its 
lifting to a certain height, calculated by the program, by rotating the drum by 180 degrees and then stopping. The height 
of the load suspension indicated on the graph is 0.78 m. 

The second curve on this graph shows the change in the lifting speed of the load over time. When the drum rotates, 
the lifting speed of the load initially increases, then gradually decreases until it stabilizes. 

The third curve represents the change in the acceleration of the load lifting over time. At the moment of lifting the 
load, the acceleration reaches a maximum value of 0.26 m/s? in 0.8 seconds, and then stabilizes. 

From the analysis of these three curves, it can be noticed that the stabilization time of the crane operation is almost constant. 

Coil Angle Curve 

Figure 5 shows the dependence of the winding rotation on time in degrees. 


Coordinate string, 6° 
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Fig. 5. Coil Angle Curve 
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Figure 5 shows the change in the angle of rotation of the drum over time. The rotation angle stabilizes when reaching 
a value of 180 degrees, after which it remains stable, meaning it repeats. 

Discussion and Conclusion. Analysis of the above graphs leads to the following conclusions. 

The mathematical model and algorithms allow for a detailed study of the motion of a bridge crane at all stages. 
Adjusting the winch operation and improving the metal structure of the crane have many positive aspects, including 
reducing dynamic displacements in the metal frame of the crane and transmission elements (such as clutch, gearbox, 
motor, and pulleys), as well as in the hoisting ropes. This reduction in dynamic loads leads to a decrease in rapid wear of 
these elements. Additionally, cost savings on maintenance and the development of an optimal metal structure design are 
achieved, as well as an increase in the crane's service life. This is evident from changes in the center of gravity of the 
load, speed, and acceleration during its lifting over time. It has been established that the height of the load during lifting, 
which is directly related to the length of the hoisting ropes, as well as the mass of the lifted load and the design of the 
main bridge, including its shape and dimensions, significantly influence the dynamic loads experienced by the structure. 

The mathematical model offers a new approach to improving the efficiency and safety of industrial operations by 
providing a more precise understanding and accounting for dynamic loads during the start-up of a bridge crane. 
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Abstract 

Introduction. Safety of navigation and development of underwater mineral deposits require the accurate detection of 
various underwater objects. The literature discusses the issues of tracking their motion and trajectory. Sonar methods are 
proposed to maintain high accuracy of underwater object positioning. High accuracy of the bearing of stereo sensors with 
an ultrashort base is noted. However, this equipment is sensitive to the sampling rate of the signals, which causes 
“sampling noise”. There are no publicly available publications dedicated to the solution to this problem. The presented 
study is designed to fill this gap. This work is aimed to study the possibility of obtaining data clarifying information about 
the bearing of underwater objects through using the phase information about echoed probing signals and an additional 
procedure for resampling the source data. 

Materials and Methods. The location of the object was determined using the experimental complex for studying 
hydroacoustic sensors created by V.A. Shirokov and V.N. Milich at the Udmurt Federal Research Center, the Ural Branch 
of the Russian Academy of Sciences. A stereo sensor with a small base (30 mm) was used compared to the distance to 
the object (~800-900 mm). Digital filtering methods and mathematical apparatus of correlation analysis of return 
hydroacoustic signals obtained by the phase method were used for data processing. 

Results. The results of comparing two methods for determining the bearing on an object are presented: by the difference 
in the time of arrival of the pulse-leading edges and by the maximum of the cross-correlation function (CCF). The change 
in bearing as the object moves, is graphically shown. The use of the leading edge of the signal caused small outliers of 
values along the entire bearing curve (less than 0.12 rad). At the maximum CCF, emissions were recorded only in some 
areas, but they were quite significant (about 0.17 rad). It showed how to select points corresponding to a smoother and 
more valid object trajectory, and how to work with erroneous points. The presented method of error correction can be 
implemented programmatically. With a quasi-harmonious signal, rare measurements of the original signal are interpolated 
by frequent calculated values. Thanks to this virtual increase in the sampling rate (oversampling), intermediate indicators 
can be recorded in the digitized source data. Interpolation of the signal values by a cubic spline allowed us to obtain 20 
points for | period of the signal instead of 5 points in the original version. In this case, the trajectory formed with the 
maximum CCF is more correct. 

Discussion and Conclusion. The direction-finding problem can be solved with the accuracy required for practical 
application. Taking into account the factor of smoothness and continuity of the object's trajectory makes it possible to 
qualitatively correct the selection of the maximum of the cross-correlation function of the stereo sensor signals. The 
proposed methods have great potential for the development of underwater vision systems. 
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AHHOTalna 

Beedenue. be3omacnoctb cyfOxoOACcTBa UH paspaOoTOK NOABOAHbIX MECTOPOXKACHHM MOMe3HbIX HCKOMAeMbIX TpeOyIOT 
TOYHOTO OOHapy2KCHHA Pa3IM4HBIX MOABOAHBIX OObeKTOB. B sHTepaType paccMaTpHBaIOTCA BOMPOCHI OTCICKUBAHHA UX 
TlepeMeLIeHHH HW TpaeKTOpHH ABwKeHUA. IIpezNararoTca MeTOAbI THAPOMOKANHU, OOecHedHBAaIOMIHe BbICOKYIO TOUHOCTb 
TO3HIMOHMPOBaHHA MOABOAHBIX OObEKTOB. OTMedeHa BbICOKaA TOUHOCTH MeseHra CTepeCOAaTIHKOB C YIbTpakOpOTKON 
6a30H. OgHako Takoe oOopyfOBaHve 4YBCTBHTeIbHO K YacTOTe JMCKpeTH3allMH CHTHaJIOB, YTO BbI3bIBAeCT «LYM 
TucKpeTu3aunM». B oTKpBITOM JocTyne HeT MyOMKalH, MOCBAILCHHBIX PelueHHIO 9TOM mpoOsempI. IpezctaBreHHoe 
YcCIeqOBaHve IpH3BaHO BOCHOJHUTb JaHHbI mpober. Lien padoTsl — u3y4eHHe BO3MOXKHOCTH NOJYYeHHA TaHHBIX, 
YTOUHAIOWMX HHPOpMalvIo O TMeweHTe MOABOAHIX OOBEKTOB 3a C4eT HCHONb30BaHHA (ha30BOH HMHPopMallHuH 
OTPaKeHHBIX 30HAMPYIOWMX CHTHAIOB UH JOMOJHUTeIBHOM NpOueAypbl MepeAMCKpeTH3alH HCXOAHBIX JaHHBIX. 
Mamepuanoi u memoooi. Mectononoxenue OObeKTAa OMpe eA C MOMOLI[bIO IKCICPHMCHTAIBHOTO KOMIWIeKca VIA 
YccueqOBaHHA THApoakycTwyeckux aTankos, co3qaHHoro B.A. LWupoxosppim u B.H. Munna B YyamMyprcKoM 
(peyepaIbHOM HCCIeqOBaTeIbCKOM WeHTpe YpasbcKoro oTyemeHua Poccniickoi akagemuu HayK. VUcnomp30Baiu 
cTepeoyqaTuHk c Maso Oa30% (30 MM) lo cpaBHeHMrO c paccTosHuem JO OObeKTa (~800—900 mm). JI oOpaboTKH 
JaHHbIX TIPHMeHATH MeTOAbI UMdppowot usbTpalwH WM MaTeMaTWYeCKHM alllapaT KOppesIAWHOHHOTO aHasIn3a 
OTpaKeHHBIX TH poakKyCTH4YeECKHX CHTHAJIOB, HOJYYCHHBIX (a3sOBbIM MeTOOM. 

Pe3yivmamot ucciedoeanua. UpeyctaBseHbl ATOM COMOCTaBJICHHA ABYX ClOcoOoB onpeyeseHHA WeeHTa Ha OOBEKT: 
10 pa3HOCTH BPpeMeHH MIpHxoa MWepeqHHUxX PPOHTOB HMITyJIbCOB HU IO MaKCHMYMY KPOCC-KOppeJIAWMOHHOM PyHKUHH 
(KK®). [papuyeckn moka3aHo u3MeHeHHe elena pu ABMKeHHH OObeKTa. Ucnomb30BaHve MepegqHero (pouta 
curHasia O6ycOBHIO HeOObIUIMe BEIOPOCHI 3HaYeHHH BAO BCel KpuBO Mesenra (MeHee 0,12 pan). [pu mMakcumyme 
KK® ppiOpocs (PUKCHPOBAJINCh JIMIb B HEKOTOPBIX OOACTAX, HO OBI JOBOJIBHO 3HA4HTeJIbBHBIMH (OKONIO 0,17 pay). 
Tloka3aHo, Kak BbIOpaTb TOUKH, COOTBETCTBYIOMNMe Oosee TiIaqKOl UM BaIMAHOM TpaeKTOpHH OOBEKTA, H KaK paOoTaTb Cc 
omMOouHBIMH TouKaMH. IIpeycTaBIeHHbI MeTOA ycTpaHeHua OWIMOKH MO%KHO peasM30BaTb Mporpammuo. IIpu 
KBa3HTapMOHHYHOM CHTHale peqKHe U3MepeHHA MCXOMHOTO CHTHasla HHTepMONMpy!IOTCA YaCTbIMH BbIMMCJICHHbIMU 
3HaueHHaAMH. baaroflapa TaKOMy BHPTyaJIbHOMY YBeJIMYeHHIO YaCTOTh! JHCKpeTH3alHN (MepeqUcKpeTH3al[HH) MO%KHO 
(PHKCHpOBaTb NpOMexKyTOYHIe WoKazaTeIM B OWMPpOBAaHHBIX MCXOAHBIX WaHHbIx. VaTeprouauHa 3HavyeHHi CurHasa 
KyOuyecKuM CraMHOM NosBoOsHa WOy4YnNTb 20 ToueK Ha | MepHoy curHasia BMecTO 5 TOYeK B HCXOJHOM BapuaHTe. B 
93TOM CJIyyae OoIee KOppeKTHa TpaeKTOpHaA, ccbOopMHpoBaHHasd C MaKCHMyMOM KK®. 

O@cystcoenue u 3akmo4uenue. 3afa4y MesneHTaljMH MOKHO PelIMTb C TOUHOCTHIO, HEOOXOAMMOM [JIA IpakTH4YecKoro 
lipHMeHenua. Yuet tbakTopa ryiaKOCTH MH HeIpepbIBHOCTH TpaeKTOPHM ABYWKeCHHA OOBEKTAa MO3BOJIAeT KaYeCCTBCHHO 
KOPpeKTHpoBaTb BbIOOp MaKCHMyMa KPOCC-KoppeJIAUMOHHOM PyHKUHH cHrHaloB cTepeogaTuuka. IpeyoxKeHHbIe 
MeTObI OOaarOT OOMSIIMM MOTECHIMAJIOM JIA paspaOoTKH CHCTeM MOABOAHOTO BUAeCHHA. 


KsrroyesBble c10Ba: OllpeseseHve MeCTONOIOKeHHA OOBEKTA B TUApocpesle, MeseHralua (Pa3z0BbIM MeTOOM, TepesHHe 
(PpOHTHI HMITYJIBCOB, KpOCcc-KOppeIAMHOHHad PYHKUMA, WYM WHcKpeTu3altHu 


Baarogapuoctu. ABTopbl BbIpakaloT OaroqapHOCTb OTNMYHbIM HHKeHepaM H.H. 3Bepepy u A.C. Tamy3uuy 3a 
TIOMOLb B CO3{aHHU TeXHHYecKuXx cpezcTB. UccreqoBaHHA, BbIMOJHeHHbIe OarofapsA Ux padoTe, cTaIH OCHOBOM JIA 
HallucaHHsd 3TOM CTaTbH. 
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OObeKT C HCHONb30BaHHeM (ha30BOH HMHPopMalHun AMpdpepenyuatbHoro cTrepeogaTanKa. Advanced Engineering 
Research (Rostov-on-Don). 2024;24(2):198—206. https://doi.org/10.23947/2687-1653-2024-24-2-198-206 
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Introduction. Safety of navigation and work in underwater mineral deposits require high-quality detection of 
underwater objects and tracking their movements [1]. Hydroacoustic sensors are used in underwater surveillance 
systems. They capture the signal reflected from the object and provide for the calculation of its location using the 
trilateration method [2]. In this process, each sensor provides information about the time interval of the probing signal 
reflected from the object. In the case of using several sensors [3], it becomes possible to solve the problem of spatial 
resection and determine the coordinates of the observed object [4]. To increase the accuracy of measurements, it seems 
promising to use stereo sensors with an ultrashort base as a receiver [5], which allow obtaining phase information [6] 
and determining the bearing on an object [7]. The possibilities of using sonar data on remote underwater targets, as 
well as on the bearing of underwater objects [8] through the use of phase [9] or frequency [10] information of reflected 
probing signals have been studied. 

At the same time, it is known that the signal sampling procedure causes errors in determining the parameters of the 
trajectory of underwater objects [11]. This applies to both distance measurements and bearing determination procedures. 
Digitization of the analog sensor signal, which is required for further digital processing, introduces the so-called 
“sampling error’. The error caused by the signal amplitude quantization is estimated by the number of quantization levels 
of the analog-to-digital converter. The error caused by time sampling is proportional to the range of the quantization time 
interval and the stress rate. Therefore, the sampling frequency of the signal is sought to be as high as possible above the 
upper frequency of the signal itself. However, increasing the sampling rate in a number of cases is limited by the 
capabilities of recording equipment: processor, analog-to-digital converter, data transmission channels, data storage 
devices. The limited possibilities of sampling in time when performing trajectory measurements do not allow for accurate 
recording of the extreme values of the trajectory (range and bearing). As a result, the accuracy of measurement results 
decreases, and outliers appear on the fixed trajectories. To eliminate this disadvantage, it is necessary to study the potential 
of oversampling the digitized sonar signal. There are no publicly available publications dedicated to solving this problem. 
The materials of this article fill in the existing gap. 

The objective of the presented study is to evaluate the possibilities of determining the bearing on an object using phase 
information, and to develop a phase method of bearing using the resampling of a digitized sonar signal. 

Materials and Methods. An algorithm for determining coordinates using the direction-finding phase method is 
considered. The results of its testing for an object moving along a circular trajectory in an experimental pool are presented. 

Setting up an experiment. During the experiments, an experimental assembly created by V.A. Shirokov and 
V.N. Milich at the Udmurt Federal Research Center of the Ural Branch of the RAS was used to determine the 
location of an object in the hydroenvironment with the direction finding by the phase method [12]. This is a 
laboratory measuring complex with a linear aqua medium in the form of an extended cylindrical tank (hydro 
waveguide) and a basin equipped with a system for generating test sonar signals. The installation elements used in 
the experiments are described in detail below. 

1. Hydroacoustic signal emitter with known coordinates S (Fig. 1). An amplitude-modulated signal is used in the work [13]. 


v 


Fig. 1. Location of the receiving sensors (stereo sensor with two receivers st] and st2) and emitter (S) in the plane of experimental 
tank (xy). Here, d — stereo sensor base; D(O; st) — distance between the stereo sensor and the object; 
o. — bearing on the object; plate (x 'y’) is formed by the plane of the stereo sensor and normal z to it 
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2. Reflecting object. In the experiment, an extended cylindrical object is used as a test object — a copper wire whose 
diameter (0.15 mm) is significantly less than the acoustic wave length (1.5 mm). In the presented work, the problem is 
solved on a plane; therefore, for the compilation of computational procedures for calculating coordinates, this object can 
be considered a point. 

3. Stereo sensor that converts a sonar signal into an electrical one. The coordinates of each of the two receivers of the 
stereo sensor are known (see sf1, s¢2 in Fig. 1). 

4. Hardware-software complex that performs amplification, digitization and processing of sensor signals. 

The measurement result is a recording file of the received two-channel signal. 

Algorithm for processing measurement results. Using a stereo sensor with a small base (30 mm) compared to the 
distance to the object (~ 800-900 mm) allows us to apply simplified formulas. 

Now then, we calculate the coordinates of the object under study in the system x'y', formed by the plane of the stereo 
sensor and the normal to it. To do this, we use the algorithm that includes five stages. 

1. Preprocessing: 

— removal of the blind zone (the first N counts); 

— signal filtering. 

High-frequency filtering is used in the work. 

2. Selecting informative fragments imp1, imp2 (in two channels) containing pulses reflected from the object under 
study. To do this, thresholding is performed according to amplitude A of the signal. The threshold is assumed to be equal 
to pxmax(A) (p = 0.5). 

This processing allows for determining the leading edge of the pulse (BeginIndex), reflected from the object. Sections 
[BeginIndex-0.5xLenSignal; BeginIndex+LenSignal] are cut from two signal channels. Here, LenSignal — length of the 
input signal. As a result, we get two signals imp1 and imp2 with a length of 1.5xLenSignal. Each of them contains the 
pulse reflected from the object. 

3. Calculation of distance D(O; st) between the object and the stereo sensor. 

— Determination of the distance traveled by the pulse from the emitter to the object and from the object to the reception 
sensor, using formula D(S,O,st)=BeginIndexxdtxC. Here, C = 1475 m/s — acoustic wave velocity, dt = 0.2x10-6 s — 
sampling interval. 

— Calculation of value 
D(S,O,st)— D(S,st) 


D(O,st) = 5 (1) 
4. Determination of bearing o on object [14]: 
a=aresin{ FAA) (2) 


where 4 — wavelength; F — signal frequency; dt — sampling interval; m — difference in the arrival of the reflected 
pulse to two stereo sensors (in counts); d — base of the stereo sensor (in the presented study 1 = 1.5 mm, F = 1 MHz, 
dt = 0.2 us, d= 30 mm). 

5. Calculation of the coordinates of the analyzed object in the system formed by the plane of the stereo sensor and its 
normal (Fig. 1), according to the formulas: x’ = D(O, st) xcos(a); y’ = D(O, st) <sin(a). 

Determining the bearing on an object. Assume the condition of the smallness of the stereo base compared to the 
distance from the sensor to the object. In this case, the stereo sensor allows us to determine the bearing on the object using 
information about the phase difference of the received pulses previously reflected from the object, as well as the difference 
in the time of arrival of the leading edges of the pulses received by the stereo sensor (Fig. 2). 


A, Va 


Fig. 2. Determination of difference in the moments of arrival of stereo signals ds, which are indicated by solid imp1(t) and 
dotted imp2(t) lines. Th1, Th2 — thresholds for the first and second stereo pair signals, respectively; dt — sampling interval; 
A — signal value in volts; t — time 
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To implement this approach, the following procedures are performed. 

1. Identification of the moments of arrival of the leading edges of stereo signals using thresholding [15] with threshold 
pxmax(A). In this paper, coefficient p is assumed to be 0.5 [16]. 

2. Calculation of difference ds (number of intervals dt) between the leading edges of the pulses registered in two 
receivers of the stereo sensor. 

3. Calculation of bearing a by formula (1). We take value d for m. 

With this method of calculating the phase shift, the leading edge of the signal is the main indicator of the pulse reflected 
from the object. The internal structure of the signal is not taken into account, which may be crucial for the correct 
determination of the bearing on the object. We take into account the similarity of the received signals and consider another 
approach to calculating the bearing — based on the cross-correlation function (CCF) of stereo signals. It is logical to 
determine the phase shift by the CCF maximum. The elements of the bearing calculation algorithm using CCF are 
described below. 

1. Calculation of the CCF pulses arrived at two receivers of one stereo sensor: r= xcorr(imp1, imp2). The value of 
function xcorr at the shift of m of the second signal relative to the first is calculated as follows: 

N=m-1 
imp\ nim X imp2,, m= 0, 


XCOMT yp, imp2 (1) Fe 70. 


XCOMmy2,imp1(—m), m<0. 

Here, N — length imp1 and imp2, —N<m<N. 

2. Determination of the maximum CCF and the corresponding shift dc (in counts). 

3. Calculation of bearing a from formula (1) taking into account that value dc is taken as m. 

Research Results 

Testing the bearing determination methods. The presented methods for determining the bearing on an object (by 
the difference in the time of arrival of the pulse-leading edges and by the CCF maximum) were tested for an object moving 
along a circular trajectory in the experimental pool. A harmonic signal with a duration of 7 periods and a frequency of | 
MHz was used. Digitized stereo signals were recorded at 256 points of the object's trajectory. The digitization frequency 
was 5 MHz. 

Figure 3 shows the change in bearing values as the object moves around the circle for the two proposed approaches 
(difference of the leading edges, CCF maximum). 


a, rad 
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Fig. 3. Deviation of bearing a on the object making a circular motion. Indicators were obtained using two approaches 
(red marker — difference in the time of arrival of the leading edges of stereo signals, blue marker — CCF maximum). 
MN — measurement number, o — bearing on the object 
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When using the signal-leading edge (Fig. 3, red marker), small outliers of values are observed almost along the entire 
curve of the bearing on the object. When using the CCF maximum (Fig. 3, blue marker), only in some areas there are 
outliers of bearing values, but they are quite large. Perhaps this is due to the peculiarities of the object's movement in the 
area where the probing and reflected signals cross the zone of turbulence caused by the object's movement. When finding 
the phase shift of stereo signals using the CCF maximum, additional errors occur, which is then reflected in the calculated 
trajectory of the object. Also, outliers in the calculated bearing values may be due to an insufficiently high sampling rate. 
The values of these errors are significantly higher than the emission values observed for the method using the pulse- 
leading edge. 

Figure 4 a presents the trajectories of the object, calculated using the global CCF maximum and its neighboring 
local maxima. 
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Fig. 4. Analysis of the trajectories of the circular motion of the object: 
a — trajectories calculated using the global CCF maximum (black marker), 
local nearest right CCF maximum (orange marker), local nearest left CCF maximum (green marker); 
b — example of CCF, whose global maximum corresponds to the correct value of the coordinates of the object; 
c — example of CCF, in which the left local CCF maximum corresponds to the correct value of the coordinates of the object 


Using aprior information about the smoothness of the movement of the experimental object, it is possible to select 
points that correspond to a smoother and more valid trajectory. This approach is used in other works [17]. In Figure 4 a, 
at x’e[800;870], vy’ >0, the correct values are on the curve for the trajectory calculated from the left local CCF maximum. 
There are also erroneous points. If they are replaced by the points of the trajectory obtained as a result of using the right 
local CCF maximum, the correct value of the trajectory point can be calculated (x ’~800, y’<0 in Fig. 4 a). This graphical 
method of error correction can be implemented programmatically. 

Thus, at some points of the trajectory, the correct phase shift of the stereo signals required for calculating the bearing 
may correspond not to the global CCF maximum, but to one of the neighboring local extremes (Fig. 4 5, c). This shift of 
the CCF maxima may be due to an insufficiently high sampling rate of the signals. This also explains the stepwise nature 
of the resulting trajectory. 

Taking into account the quasi-harmonicity of the signal, it is proposed to use interpolation of rare measurements of 
the original signal with frequent calculated values. This virtual increase in the sampling rate (oversampling) makes it 
possible to fix intermediate values in the digitized source data. Interpolation of the signal values by a cubic spline allowed 
us to obtain 20 points for | period of the signal instead of 5 points in the original version. Figure 5 shows the calculated 
trajectories of the object, formed with preliminary resampling of the signals. 
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Fig. 5. Calculated object motion trajectories obtained for stereo signals with a virtually increased sampling rate: 
blue marker — difference in the time of arrival of the signal-leading edges is used; 
red marker — CCF maximum is used 


The red line is the object’s trajectory, obtained using the CCF maximum to calculate the bearing from pre-interpolated 
data. It seems to be more correct compared to a similar trajectory obtained using the difference in the leading edges of 
reflected pulses (blue line). 

Discussion and Conclusion. Analysis of the features of determining the bearing on an underwater object using the 
mutual phase information of the differential stereo sensor signals makes it possible to draw the following conclusions. 

— The difference in the arrival time of the pulse-leading edges recorded by the stereo sensor allows us to obtain bearing 
values indicating the direction of the object. The quality of the result depends on the accepted threshold value (used in 
determining the pulse-leading edge) and the variability of the amplitudes of the received signals. 

— The application of the approach using the maximum of the cross-correlation function at an insufficient sampling 
rate of the initial signals leads to significant outliers in the calculated trajectory (Fig. 3). 

— Resampling of signals through interpolation of the original quasi-harmonic signals using the CCF maximum for the 
calculation of the bearing and coordinates of the object provides obtaining a smooth trajectory (Fig. 5), corresponding to 
the smooth movement of the object in a circle. This approach requires significantly higher computational costs compared 
to the others discussed in this article. This may affect the speed of real-time signal processing. 

Thus, the investigated possibilities of refining the phase direction finding method (resampling of digitized signals and 
using the CCF maximum) make it possible to effectively assess the trajectory of an underwater object. 

The conditions for obtaining a smooth correct trajectory of the tracked object are described. The oversampling of 
interpolated quasi-harmonic stereo signals digitized with insufficient sampling rate serves this goal. At the same time, the 
bearing calculation method is used with the maximum of the cross-correlation function of the stereo sensor signals. 

The data in Figure 5 confirm that the direction-finding problem can be solved with the accuracy required for practical 
application. Taking into account the factor of smoothness and continuity of the object's trajectory allows for the unique 
adjustment of the selection of the maximum of the cross-correlation function of the stereo sensor signals (Fig. 4). 

The proposed methods are of great importance for the development of underwater vision systems. 
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